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ABSTRACT 


Until  the  late  sixties  and  seventies,  the  main  approaches  to 
treating  fluid  dynamics  in  the  turbulent  regime  were  based  on 
statistical  techniques  and  semiempirical  models.  The  advancements  in 
applied  and  nonlinear  mathematics  initiated  research  in  the  stability 
properties  of  equations  governing  fluid  motion.  This  allowed  the 
modelling  of  the  breakdown  of  certain  flow  structures  into  the  chaotic 
state.  Physical  models  describing  possible  mechanisms  of  turbulence 
in  connection  with  the  aforementioned  stability  properties  of  the 
governing  equations  are  still  rare.  The  present  work  describes  such 
a  model  and  its  machine  simulation. 

The  mechanism  for  the  generation  of  turbulence  assumes  a 
sequence  of  three  Hopf-bifurcations  on  the  Navier-Stokes  equations. 
This  leads  to  a  3-torus  as  invariant  manifold  in  state  space.  The 
vector  fields  on  this  3-torus  are  structurally  unstable  and  chaotic 
behaviour  is  triggered  by  combinations  of  small  random  perturbations. 
To  simulate  this  mechanism,  two  Navier-Stokes  type  equations  for  the 
evoluation  of  time  dependent  velocity  perturbations  are  derived  and 
discretized  in  three  dimensions  by  using  finite  elements  and  the 
Galerkin  method.  The  resulting  systems  of  ordinary  differential 
equations  are  then  tested  by  a  FORTRAN  program  for  Hopf  bifurcations 
and  solution  curves  in  state  space  were  obtained  numerically  on  the 
computer.  Small  random  perturbations  are  imposed  on  these  solution 
curves  by  using  random  number  generators  and  a  spectral  analysis  for 
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energy  densities  using  Fast  Fourier  transform  under  use  of  the  Taylor 
hypothesis  is  performed  on  several  of  these  solutions  with  and  without 
random  perturbation  in  order  to  obtain  energy  spectra  which  permit 
comparison  with  the  inertial  subrange  law  for  homogeneous  and  isotropic 
turbulence. 
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INTRODUCTION 


Perhaps  one  of  the  main  reasons  why  the  phenomenon  of  turbulence 
lacks  any  satisfactory  description  or  model  is  the  fact  that  no  completely 
satisfactory  definition  of  turbulence  seems  to  be  possible.  One  of  the 
simplest  and  straightforward  definitions  might  be  'a  state  of  continuous 
instability*  [Tritton,  [1]].  This  again  can  be  understood  only  intui¬ 
tively,  because  the  concepts  of  continuity  and  instabil ity  are  almost 
contradictory;  continuous  processes  take  place  over  a  considerable 
length  of  time,  while  the  phenomenon  of  instability  occurs  as  a  sudden 
breakdown  of  a  state  and  the  following  transition  into  another  one. 

Looking  at  turbulence,  one  cannot  observe  different  states  at  different 
points  in  time,  if  the  governing  external  parameters,  such  as 
temperature,  temperature  gradient,  mean  shear,  or  viscosity  are  not 
changed.  Moreover,  it  is  not  certain  whether  there  exists  only  one  main 
mechanism  in  nature,  which  governs  all  the  different  types  of  observed 
turbulence . 

So  far,  only  methods  of  statistical  analysis,  which  lead  to  the 
various  techniques  of  finite  closure  of  moment  equations  derived  from 
the  Navier-Stokes  equations*  have  given  results  which  are  comparable  with 
experiments.  One  of  the  most  outstanding  results  in  this  respect  is  the 
inertial  subrange  law  first  derived  by  Kolmogorov  using  dimensional 
analysis  on  a  rather  heuristic  model  of  energy  transfer  within  isotropic 
turbulent  flow.  Using  methods  of  statistical  physics  and  finite  clo¬ 
sure  techniques  applied  to  these  moment  equations,  Heisenberg  reproduced 

• k 

They  will  be  referred  to  as  NS-equations  from  now  on. 
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this  result  which  can  be  verified  experimentally  (see  Appendix  E) .  It 
is  sometimes  also  referred  to  as  the  (-5/3) -law,  since  it  states  the 
proportionality  between  spectral  energy  density  E  and  wave  number  k 
to  the  power  of  -5/3  : 


For  an  introduction  to  this  result  see  e.g.  [1]  or  Appendix  E. 

The  application  of  methods  from  differential  topology  and 
ergodic  theory  to  the  qualitative  theory  of  differential  equations* * 
lead  to  the  development  of  concepts  such  as  structural  stability  and 
strange  attractors.,  and  thereby  to  a  better  understanding  of  their 
solution  spaces.  However,  due  to  the  difficulties  encountered,  only 

very  slow  progress  can  be  made  even  for  low-dimensional  systems  of 
ordinary  differential  equations.**  Complete  results  about 

structural  stability  exist  so  far  only  for  two-dimensional  (2D)  systems 
and  the  term  'strange  attractor'  was  coined  to  reflect  the  lack  of 
understanding  of  the  phenomenon  it  tries  to  depict.  This  term  is  used 
mainly  in  connection  with  chaotic  behaviour  of  solutions  of  differential 
equations  and  the  attracting  regions  in  solution  space  associated  there¬ 
with.  In  the  context  of  fluid  dynamics  chaotic  behaviour  is  called 
turbulence;  again  one  experiences  the  lack  of  ability  to  come  to  a 
satisfactory  definition  and  has  to  resort  to  terms  such  as  'strange 
attractor' . 


For  an  introduction,  see  e.g.:  Nemytskii  and  Stepanov:  Qualitative 
theory  of  differential  equation.  Princeton. 

*  Jc 

They  will  be  referred  to  as  ODE-systems  from  now  on. 
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Since  the  set  of  NS-equations  is  far  from  being  a  low-dimensional 
ODE-system  but  rather  a  system  of  three  nonlinear,  parabolic  partial 
differential  equations*,  its  state  space  is  infinite-dimensional. 
Moreover,  its  solution  type  can  change  drastically  with  the  boundary 
conditions  or  the  geometry  of  fluid  domain  imposed. 

As  a  consequence,  the  analysis  of  equations  of  the  NS-type  is 
being  conducted  in  two  related,  but  nevertheless  different  ways.  One 
approach,  mainly  taken  by  engineers,  deals  with  very  special  boundary 
conditions  and  tries  to  find  the  associated  flows  as  solutions  of  the 
NS-equations  which,  as  a  rule,  first  have  to  be  adapted  to  the  problem 
at  hand.  In  this  connection,  various  techniques  of  discretizing  the 
fluid  domain  can  be  used  to  obtain  numerical  solutions.  The  other  way 
is  to  approach  the  NS-equations  on  a  very  general  level  and  try  to 
analyze  the  structure  of  its  solution  space.  This  is  mainly  done  by 
mathematicians,  since  the  techniques  involved  require  extensive  use  of 
higher  analysis,  differential  topology  and  function  theory. 

One  of  the  better  known  publications  in  this  field  which  suggests 
a  mathematical  mechanism  for  the  generation  of  turbulence  and  chaos 
originating  from  solutions  of  equations  of  the  NS-type  is  the  paper  of 
Ruelle  and  Takens  [2] .  These  authors  have  shown  that  the  proposition 
of  Landau  and  Lifschitz  [3]  for  turbulent  motion  as  a  quasi-periodic 
function  in  time  is  invalid. 

The  present  study  was  largely  motivated  by  the  aforementioned 
paper  of  Ruelle  and  Takens.  Therefore,  after  an  introduction  to  the  Hopf 
bifurcation  (Chapter  1) ,  which  is  the  main  tool  for  the  construction  of 

★ 

They  will  be  referred  to  as  PDE  from  now  on. 
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the  mathematical  mechanism,  a  rather  informal  discussion  of  this  paper 
follows  (Chapter  2).  After  that,  the  model,  which  is  the  main  subject 
of  the  present  study,  and  its  mechanism  are  explained  (Chapter  3) .  This 
is  done,  first,  by  an  informal  description  aimed  at  explaining  the  basic 
ideas  in  terms  as  simple  as  possible.  Following  this  rather  heuristic 
introduction  is  a  mathematical  formulation  of  the  same  model.  The 
necessary  function  spaces  are  introduced,  which  permit  representing  of  the 
model  by  means  of  equations.  This  concludes  the  first  part  and  there¬ 
with  the  description  of  the  theoretical  basis. 

Part  2  describes  the  technical  realization  of  the  ideas 
described  in  Part  One.  First,  a  set  of  equations  of  the  NS-type  is 
developed  for  a  set  of  perturbations  on  the  fluid  velocity  superimposed 
on  each  other  (Chapter  4) .  To  obtain  numerical  solutions  for  this  set 
of  equations,  a  discretization  of  the  domain  in  spaces  governed  by  these 
equations  is  carried  out,  based  on  finite  element  techniques  and  the 
Galerkin  method  (Chapter  5) .  This  leads  to  an  approximation  of  NS-type 
PDE's  by  a  system  of  ODE's.  In  connection  with  this  discretization, 
divergence-free  basis  functions  or  finite  elements  (Chapter  6)  and  the 
use  of  Dirichlet  boundary  condition  for  the  domain  of  the  velocity 
perturbations  are  introduced  and  a  particular  mean  flow  structure  is 
discussed  (Chapter  7) . 

Part  3  describes  the  numerical  study  of  the  system  of  ODE's 
obtained  through  the  discretization  discussed  in  Part  2.  First,  the 
system  of  ODE's  is  investigated  for  its  bifurcation  properties  for 
different  profiles  of  the  underlying  mean  flow  (Chapter  8) .  The  results 
of  four  particular  mean  flows  are  selected  for  further  analysis  of  their 
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stability  characteristics  by  use  of  numerical  ODE-solvers  and  graphs 
of  their  solution  trajectories  (Chapter  9).  One  of  these  four  flows 
is  studied  under  the  superposition  of  numerically  simulated  random 
perturbations  (Chapter  10)  and  the  results  are  submitted  to  a  spectral 
analysis  to  investigate  any  realization  of  the  inertial  subrange  law 
(Chapter  11) . 

Finally,  the  summary  and  conclusions  (Chapter  12)  contain 
an  overall  interpretation  of  the  results  obtained  in  the  three 
parts.  Also,  a  number  of  technical  restrictions  and  their  effect  and, 
in  connection  therewith,  the  comparability  of  the  spectral  analysis 
with  experimentally  obtained  spectra  is  considered. 


V 


PART  1 


1.  THE  HOPF  BIFURCATION 

The  following  introduction  to  the  Hopf  bifurcation  is  rather 
informal  and  concise.  Although  Ruelle  and  Takens  apply  it  to  NS-type 
PDE’s  it  is  being  introduced  here  for  ODE’s  for  the  sake  of  simplicity. 
More  elaborate  descriptions  with  theorems  and  proofs  can  be  found  in 
various  publications,  see  e.g.,  Marsden  and  McCracken  [4]  or  Hassard, 
Kazarinoff  and  Wan  [5]  or  Hopf's  original  paper  [6], 

Consider  an  n-dimensional  ODE-system  (n  >_  2)  : 

X ' (t)  =  f(X  (t),u)  (1.1) 

For  the  real  parameter  \±  within  a  certain  interval  I,  the  system  is 
assumed  to  have  the  stationary  solution  X^s  according  to  X^(t)  =  0, 
i.e. , 

0  =  f(X  u),  U  S  I.  (1.2) 

y  ^ 


Taking  the  variational  derivative  of  (1.1)  at  X  yields  the 

yL 

linearized  system: 


X’(t)  =  A(y  )X  (t); 

M  u  M 


A(liL)  E 


3f(Xii(t),pL) 

3X  (t) 


X  (t)=X 

y  y. 


(t) 


(1.3) 


A(y^)  is  the  Jacobian  matrix  of  the  system  at  X  (t)  =  X^  (t) .  A  Hopf 

L 

bifurcation,  i.e.,  the  appearance  of  a  periodic  solution  X  (t)  which 
grows  in  amplitude  with  increasing  or  decreasing  y  out  of  the 
stationary  solution  at  y  =  ,  occurs  when  the  following  conditions 
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are  satisfied: 


(I) 


(II) 


(III) 

(IV) 


A  pair  of  complex  conjugate  eigenvalues  X(y  )  and  X(y  )  of 

c  c 

A(pc)  cross  the  imaginary  axis  at  y  =  yc  . 

All  the  remaining  eigenvalues  of  A(y  )  stay  away  from  the 

Li 

imaginary  axis  for  y  near  y^  . 

Re  X(yc)  =  0;  Im  X  CyQ)  ^  0  • 


dX(y) 

dy 


t  0  . 


y  =y 


c 


Properties  of  stability  and  criticality  of  the  periodic  solution,  i.e., 

on  what  side  of  y  X  ft)  exists  and  whether  it  is  stable  or  unstable 

c  yp 

require  a  more  detailed  analysis  of  the  nonlinear  terms  in  Equation 

(1.1). 


2.  THE  WORK  OF  RUELLE  AND  TAKENS 

In  their  paper  Ruelle  and  Takens  consider  a  sequence  of  Hopf 
bifurcations  occuring  to  equations  of  the  type  (1.1).  However,  they  do 
not  restrict  themselves  to  ODE's  and  include  PDE's  expecially  of  the  NS- 
type.  The  parameter  y  can  therefore  be,  e.g.,  mean  velocity, 
viscosity,  shear  or  quantities  like  Reynolds  number,  Rayleigh  number 
etc.  depending  on  the  type  of  application.  Moreover,  the  right-hand 
side  of  Equation  (1.1)  has  to  be  understood  as  an  operator  generating 
a  vector  field  in  state  space*  H  of  Equation  (1.1).  In  the  ODE  case 
this  'RHS-operator '  can  be  represented  as  a  combination  of  matrices  with 
not  necessarily  linear  properties  as  far  as  its  effect  on  elements  in 
H  is  concerned.  In  the  case  of  PDE’s  the  combination  consists  of 

•k 

This  space  is  also  called  phase  space  by  some  authors,  see  e.g.  [9]. 
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the  standard  differential  operators  such  as  gradient",  divergence  and 
Laplacean,  again  with  possibly  nonlinear  behaviour.  The  eigenvalues 
to  be  considered  are  now  the  ones  of  the  linearized  forms  of  these 
'RHS-operators ' . 

The  first  part  of  the  paper  contains  a  description  of  a  potential 
mechanism  for  turbulence  by  successively  advancing  pairs  of  complex- 
conjugate  eigenvalues  of  the  linearized  RHS-operators  over  the 
imaginary  axis  as  the  critical  parameter  y  increases.  As  the  first 
pair  crosses  the  imaginary  axis,  a  fixed  point  in  solution  space  (i.e. 
a  stable  stationary  solution)  exchanges  stability  with  a  family  of 
closed  Hopf  orbits  (i.e.,  periodic  solutions)  which  grow  monotonically 
with  y  .  The  crossing  of  the  second  pair  causes  this  first  orbit  to 
loose  stability  and  degenerate  into  a  trajectory  on  the  surface  of  a 
2-torus.  This  2-torus  is  an  attracting  invariant  manifold,  i.e.  the 
trajectory  of  any  solution  will  eventually  end  up  on  its  surface. 

Further  increase  of  the  parameter  y  causes  the  third  pair  of  complex- 
conjugate  eigenvalues  to  cross  the  imaginary  axis  and  the  degeneration 
of  the  2-torus  into  a  3-torus.  This  is  the  point  where  the  subject 
matter  becomes  difficult  to  understand,  since  a  3-torus  can  only  'exist' 
in  a  space  of  at  least  four  dimensions.  Problems  of  visualization  arise 
especially  with  the  geometry  of  the  trajectories  for  which  this  3-torus 
acts  as  attracting  invariant  manifold. 

The  geometry  of  the  trajectory  on  the  surface  of  the  2- torus  is 
self-evident;  it  can  be  viewed  as  a  line  wound  into  a  coil  that  encloses 
the  torus  (see  Figure  1)  .  That  means,  this  trajectory  is  determined  by  a 
vector  field  on  the  surface  of  the  2-torus  with  a  slope  or  pitch  which 
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FIGURE  1. 


2-Torus  with  Trajectory  and/or  Vector  Field  on 
with  Velocity  Components  V^ ,  V0  and  Periods 


its  Surface 

1VP2 
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is  determined  by  the  ratio  of  the  periods  p  ,  p9  arising  from  the  two 
bifurcations  (see  Figure  1) . 

Crossings  of  further  eigenvalue  pairs  over  the  imaginary  axis 
cause  further  Hopf  bifurcations  and  thereby  appearance  of  the  appropriate 
higher-dimensional  tori  as  invariant  manifolds. 

It  is  evident  that  the  mathematical  description  and  formulation 
of  the  bifurcations  following  the  first  one  is  by  far  not  as  straight¬ 
forward  as  this  very  first  one,  which  is  a  direct  application  of  the 
Hopf  bifurcation  theorom  on  the  (linearized)  RHS-operator  of  equations 
of  the  type  (1.1)*.  The  second  bifurcation  does  not  grow  out  of  a 
fixed  point  or  a  stationary  solution,  but  one  that  is  already  periodic 
in  time.  Ruelle  and  Takens  therefore  introduce  a  Poincare  map  defined 
on  a  piece  of  hypersurface  transversal  to  the  first  Hopf  orbit  (see 
Figure  2) .  The  penetration  point  of  the  first  Hopf  orbit  on  the 
Poincare  hypersurface  is  naturally  a  fixed  point  for  the  Poincare  map. 

It  might  of  course  change  position,  if  the  Hopf  orbit  changes  size  and 
shape  due  to  variations  in  the  bifurcation  parameter  p  .  However,  if 
this  first  orbit  destabilizes,  so  will  the  fixed  point  of  its  associated 
PoincarS  map.  If  this  destabilization  is  due  to  a  second  Hopf 
bifurcation,  this  fixed  point  will  degenerate  into  an  invariant  set  of 
points  sitting  on  a  closed  orbit  on  the  Poincare  hypersurface.  The 
eigenvalues  of  this  Poincare  map  are  in  first  (linear)  approximation 
the  same  as  the  Floquet  coefficients  of  the  destabilizing  first  orbit. 

ic 

In  the  case  of  PDE ' s  even  this  can  become  prohibitively  difficult, 
since  the  analysis  of  the  spectrum  of  the  RHS-operators  in  order  to 
find  functional  dependencies  between  eigenvalues  and  possible 
bifurcation  parameters  is  definitely  nontrivial. 
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FIGURE  2.  Poincare  Map  on  Plane  Transversal  to  Hopf  Orbit 


FIGURE  3.  Circle  Oscillating  within  Annulus 
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Most  of  the  technical  work  in  the  paper  consists  of  the  statements  and 
proofs  of  the  Hopf  bifurcation  theorems  for  vector  fields  (relating  to 
the  first  bifurcation)  and  diffeomorphisms  (e.g.,  Poincare  maps  used 
for  the  second  bifurcation)  and  supporting  theorems  and  lemmas. 
Bifurcations  following  the  second  one  which  lead  to  higher -dimensional 
tori  are  not  being  discussed  in  rigorous  mathematical  detail  in 
this  paper.  In  the  meantime,  however,  extensive  analytic  work  has  been 
done  in  this  direction  by  loos  and  Chenciner  [7]  and  Sell  [8] . 

Despite  the  aforementioned  difficulties  arising  with  the  complex 
structure  of  vector  fields  and  trajectories  with  higher-dimensional 
tori  as  invariant  manifolds,  there  is  no  reason  at  this  point  to  assume 
a  breakdown  of  their  deterministic  nature  and  a  transition  to  chaotic 
and  turbulent  behaviour.  -However,  in  the  last  section  of  their  paper, 
Ruelle  and  Takens  show  how  small  perturbations,  imposed  on  vector  fields 
defined  on  tori  of  dimension  higher  than  two,  can  give  rise  to  transitions 
to  vector  fields  which  are  not  Morse-Smale  (see  Appendix  1  for  defini¬ 
tion  and  explanation)  and  which  have  the  already  mentioned  strange 
attractors.  Before  this  part  of  the  paper  is  discussed,  a  few 
geometric  interpretations  of  projections  and  Poincare  maps  of  the 
3- torus  are  given. 

Advancing  the  idea  of  the  Poincare  map  from  the  2-torus  (where 

it  arranges  the  points  on  a  closed  orbit,  as  already  described)  to  a 

3-torus  one  obtains  as  invariant  manifold  for  the  Poincare  map  a  2-torus. 
Similar  considerations  can  be  used  for  projections;  the  annulus  which 

arises  from  the  projection  of  a  2-torus  into  a  plane  parallel  to  the 

plane,  which  holds  its  first  circle,*  relates  to  a  2-torus  on  which  the 

— 

See  footnote  on  next  page. 
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surface  becomes  a  shell  with  a  certain  thickness  (as  the  annulus  can  be 
understood  as  a  circle  with  a  certain  finite  'line  thickness')  if  a  3- 
torus  is  projected  into  the  three  dimensional  space  which  hold  the 
cartesian  product  of  its  first  two  circles.*  A  way  to  visualize  the 
role  of  the  dimension  lost  due  to  projection  is  to  substitute  for  it  by 
time;  in  the  case  of  the  annulus  this  leads  to  a  circle  oscillating 
radially  within  the  annulus  (see  Figure  3)  and  in  the  higher-dimensional 
case  of  the  torus  shell  one  can  imagine  a  2-torus  oscillating  within  the 
shell  (see  Figure  4).  This  approach  works  especially  well  in  cases  like 
the  present  one,  where  the  tori  are  invariant  manifolds  of  dynamical 
systems,  that  is,  time  dependent  systems,  and  it  is  evident  that  the 
oscillation  frequencies  would  have  to  be  associated  with  the  periods  of 
the  circles  or  closed  orbits,  respectively,  which  collapse  to  line 
segments  under  the  described  projections,  i.e.  the  orbits  relating  to 
the  second  or  third  bifurcation  respectively.  Obviously,  this  idea 
works  also  with  projections  into  other  planes  or  three-dimensional 
spaces  than  the  ones  holding  the  orbits.  However,  the  resulting  mani¬ 
folds  would  be  much  more  complicated  than  the  annulus  and  torus  shell, 
since  properties  of  connectedness  of  the  2-  and  3-torus  might  change 
drastically  under  such  projections.  Especially  substitution  of  time 
for  the  missing  dimension  would  lead  to  rather  complicated  motions. 

Ruelle  and  Takens  define  the  following  vector  field  z  on  the 
3-torus : 

■k 

An  n-torus  embedded  into  an  m(>  n) -dimensional  space  can  be  understood 
as  a  noncommuting  cartesian  product  of  n  circles:  Tn  =  x  ...xT 

By  the  terms  first  circle  or  first  two  circles  respectively,  or 

respectively  is  meant  here. 
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FIGURE  4.  2-Torus  Oscillating  within  Torus  Shell 


15 


0). 

oj  =  |  co 2  1  is  a  constant  vector  field  consisting  of  the  three 
w3 

angular  velocities  associated  with  the  three  Hopf-orbits,  which 
generate  the  3-torus  in  their  cartesian  product. 

/  xi  \ 

x  =  y  x  J  is  the  suspension  of  a  horseshoe  diffeomorphism  on 

the  surface  of  the  2-torus  associated  with  the  w.-  and  to^-otbits. 
That  is,  x  rearranges  the  points  of  a  disk  on  the  surface  of  the 
2-torus  into  a  horseshoe  (as  described  by  S.  Smale  in  [10];  see  also 
Figure  7)  within  the  area  originally  occupied  by  the  disk.  This 
rearrangement  is  completed  after  one  return  period  of  co^.*  In  other 
words,  the  horseshoe  diffeomorphism  is  here  simultaneously  a  Poincard 
map  associated  with  co^  •  Dividing  x  by  an  integer,  say  m,  leads 
to  x/m  .  This  means,  that  the  above-mentioned  disk  will  be  horseshoe- 
diffeomorphed  once,  after  a  sequence  of  m  return  periods  (i.e., 
Poincare  maps)  of  (in  difference  to  the  aforementioned  case,  where 

x  horseshoe-diffeomorphs  the  disk  after  one  return  period) .  For  m 
large  enough,  the  vector  field  x/m  becomes  arbitrarily  small.  This 
happens  in  z,  where  the  role  of  m  is  taken  over  by  and  q0  . 

For  q^  and  q9  large  enough,  the  difference  between  z  and  oj  can 


Clearly,  this  can  be  visualized  by  means  of  the  2-torus 
oscillating  with  period  of  co^  in  the  shell  of  the  projection  of  the 
3-torus . 


v 
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Fig.  7*  Horseshoe  d i f feomorph i sm  sketched  on  the  outer  surface  of 
the  torus  shel 1 . 
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be  made  small  enough,  to  put  z  and  oj  in  any  C  -neighbourhood  of 
each  other.  The  particular  structure  of  z,  i.e.,  the  factor 
and  the  unit  vector  in  the  third  component  are  some  of  the  requirements 
which  guarantee  the  existence  of  a  Cantor  set  in  the  set  of  non¬ 
wandering  points  of  z  despite  the  presence  of  the  part  co  .  The  other 
requirements  are:  0  <_  £  to9  <.  ;  l  <  q,/q2  <  ^  » 

0  <  w1/a)3  =  p1/q1  <  1  ;  0  <  w2/w3  =  P2/q2  <  1  ;  P1»P2»q1*q2  are 

integers  and  p^,q2,q^,p2  have  no  common  divisor.  Moreover,  the  afore- 

2 

mentioned  disk  on  T  does  not  occupy  any  points  outside  the  center 
square  which  arises,  when  the  unit  square,  defining  the  torus  by 
modulo  conditions,*  is  partitioned  into  nine  equal  squares.  These 
modulo  conditions  also  guarantee  the  formation  of  a  fi-set  which 
contains  a  set  that  is  locally  the  product  of  a  line  interval  and  a 
Cantor-set.  Since  the  vector  field  of  a  Morse-Smale  system  has  a  non¬ 
wandering  set  which  consists  of  a  finite  number  of  fixed  points  and 
periodic  orbits  only  and  since  in  the  present  system  the  Q-set  is  a 
subset  of  the  non-wandering  set,  the  present  system  cannot  be  Morse- 
Smale.  Besides,  the  fi-set  does  not  qualify  as  a  manifold  because  of  the 
above-explained  local  product  structure  and  therefore  qualifies  as  a 
strange  attractor. 

In  the  case  of  a  torus  with  more  than  three  dimensions  the 
horseshoe  diffeomorphism  is  replaced  by  a  diffeomorphism  which  double¬ 
loops  the  solid  2-torus  and  maps  it  into  the  volume  originally  occupied 

A  circle  T  can  be  identified  with  the  unit  interval  I  by  the 
mapping:  T  ->  IR  ITL  =  I  or:  t  =  r(mod  1)  ;  t  e  T;  re  IR  .  Since 

2  2 
T  =  T  x  T,  it  can  be  identified  with  I  =1x1. 
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by  it.  This  leads  to  an  ft-set  which  is  locally  the  product  of  a  Cantor 
set  and  a  piece  of  2-dimensional  manifold.  Again  this  is  the  reason 
why  it  fails  to  be  a  manifold  and  is  therefore  called  a  strange 
attractor. 

All  these  considerations  work  equally  well  if  the  vector  fields 
under  consideration  on  the  higher-dimensional  tori  are  replaced  by  ones 
sufficiently  close  to  them  (e.g.,  replace  z  by  z'  with  |z-z'| 
sufficiently  small  in  the  case  of  the  3-torus) .  The  described  cases 
are  therefore  definitely  generic  (see  Appendix  A) . 

3.  THE  MODEL  AND  ITS  MECHANISM 
Informal  Description 

Looking  at  the  infinite-dimensional  state  space  of  the  NS- 
equation  defined  over  a  certain  domain  in  fluid  space,  one  realizes 
that  it  is  spanned  by  an  infinite  number  of  orthogonal  3-dimensional 
subspaces  which  relate  to  the  infinite  number  of  points  in  the  fluid 
space  domain  with  their  attached  velocity  vectors  decomposed  into  their 
three  components.  Let  us  define  a  vector  field  over  the  NS-solution 
space  determined  by  the  accelerations  in  fluid  space  according  to  the 
NS-equations .  The  projection  of  this  vector  field  into  one  of  these 
3-dimensional  subspaces  determines  the  motion  of  the  velocity  vector  at 
the  point  in  fluid  space  associated  with  this  3-dimensional  subspace. 

In  case  the  system  has  undergone  two  stable  Hopf  bifurcations  its 
invariant  manifold  will  be  the  already  described  2-torus  and  its 
projection  into  one  of  these  3-dimensional  subspaces  will  be  diffeo- 
morphic  to  its  original  in  the  infinite-dimensional  state  space 
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(except  for  nongeneric  cases  where  it  might  appear  as  some  simply 
connected  2-dimensional  manifold,  like  the  circle  might  appear  as  a 
line  element  if  the  direction  of  projection  lies  within  its  plane)  and 
evidently  the  same  goes  for  the  vector  field  defined  on  its  surface. 

The  tip  of  the  velocity  vector  associated  with  this  subspace  will 
therefore  be  moving  on  a  trajectory  embedded  on  the  surface  of  this 
projected  torus.  The  trajectory  itself  will  be  winding  around  the  torus 
with  a  pitch  that  is  determined  by  the  ratio  of  the  two  frequencies 
(or  periods)  stemming  from  the  two  bifurcations  (see  Figure  1  ) . 

Advancing  this  projection  idea  to  the  case  of  a  sequence  of 
three  bifurcations  and  the  resulting  3-torus  in  state  space,  the 
properties  of  the  projection  into  a  3-dimensional  subspace  change 
drastically.  First  of  all,  the  projection  of  the  3-torus  from  a  space 
with  more  than  three  dimensions  into  a  space  with  three  dimensions 
cannot  be. diffeomorphic  to  its  original  (i.e.,  the  3-torus).  Recalling 
the  discussion  of  the  geometric  interpretation  of  Poincare  maps  in 
Chapter  2,  it  is  evident  that  for  the  projection  along  the  axis  of  the 
primary  orbit  into  a  3-dimensional  subspace  orthogonal  to  it  will  be 
the  already  described  torus  shell.  If  the  axis  of  the  primary  orbit 
is  not  orthogonal  to  the  3-dimensional  projection  space,  the  properties 
of  connectedness  might  change  as  well.  For  instance,  the  hollow  space 
enclosed  by  the  torus  shell  might  vanish  locally  due  to  overlapping  of 
walls.  (Again  this  is  in  analogy  to  the  lower-dimensional  case  where 
the  opening  in  the  middle  of  the  annulus  vanishes  due  to  overlapping 
effects  if  the  original  2-torus  is  being  projected  "sideways"  .) 

The  structure  of  the  vector  field  within  the  projection  of  the 
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3-torus  can  easily  be  derived  by  using  the  already  discussed  method 
of  replacing  the  dimension  lost  under  projection  by  time.  This  leads 
to  a  2-torus  oscillating  within  the  torus  shell. 

The  vector  field  on  the  surface  of  this  oscillating  2-torus  is 
evidently  the  one  already  introduced  (see  Figure  1)  and  in 
addition  to  it,  there  is  now  a  third  component,  perpendicular  to  the 
surface  of  the  oscillating  2-torus.  This  component  can  be  either  zero, 
if  the  2-torus  is  in  its  maximum  or  minimum  amplitude  position,  or 
directed  outward  (2-torus  expanding)  or  inward  (2-torus  contracting) 

(see  Figure  8) .  One  realizes  therefore  two  possible  vector  fields 
within  the  torus  shell,  according  to  the  two  possible  directions  of  this 
third  vector  component.  This  situation  can  degenerate  into  more  than 
two  different  vector  fields  if  a  projection  with  locally  overlapping 
walls  is  considered  (see  Figure  9) .  The  motion  of  the  velocity  vector 
of  a  particular  point  in  fluid  space  after  three  bifurcations  is 
therefore  considerably  more  complicated  than  after  only  two  bifurcations. 

So  far,  the  motions  resulting  from  the  vector  fields  are  still 
regularly  behaved  irrespective  of  whether  they  relate  to  two  or  three 
bifurcations.  Therefore,  a  mechanism  has  to  be  introduced  which  induces 
chaotic  behaviour. 

Looking  again  at  the  case  of  two  bifurcations  first,  one 
realizes  that  every  one  of  the  projections  into  the  3-dimensional 
subspaces  is  the  already  discussed  2-torus  with  a  local  vector  field 
structure  as  shown  in  Figure  10.  There,  a  surface  element  of  this 
2-torus  with  the  embedded  and  surrounding  vector  field  is  shown. 

Suppose  now  that  a  small  random  perturbation,  which  is  short  in  time 
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FIGURE  8. 


Vector  Field  on  Oscillating  2-Torus  with  Normal  Components 
for  Contraction  and  Expansion 
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FIGURE  9.  Vector  Fields  on  3-Torus  Section  with  Overlapping  Walls 
Region 


FIGURE  10.  Local  Vector  Field  Structure  of  2-Torus  in  3-Dimensional 


Subspace 
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as  well,  becomes  effective  at  an  individual  point  in  fluid  space.  This 
will  uncouple  the  subsystem  defined  by  this  point  from  the  rest  of  the 
system  and  lead  to  a  displacement  of  it  which  has  two  components:  one 
perpendicular  to  the  toral  surface  and  the  other  one  within  it.  After 
the  perturbation  has  relaxed  to  zero,  the  subsystem  will  recouple  to 
the  rest  system,  and  the  perpendicular  component  will  vanish,  since  the 
surrounding  vector  field  will  draw  the  system  back  into  the  stable  toral 
surface.  The  displacement  component  within  the  surface,  however,  will 
not  vanish  and  the  subsystem  will  follow  the  vector  field  on  the 
toral  surface  again,  but  slightly  displaced  sideways  w.r.t.  its 
trajectory  before  the  perturbation.  Since  random  perturbations  give 
rise  to  displacements  with  directions  equally  distributed  within  the 
toral  surface,  the  mean  trajectory  of  the  perturbed  subsystem  (averaged 
over  a  time  containing  several  random  perturbations)  has  to  be  identical 
to  the  trajectory  of  the  unperturbed  subsystem.  This  mechanism  holds 
for  all  points  in  fluid  space  and  their  state  subspaces.  The  fluid 
motion  after  two  bifurcations  must  therefore  still  be  regularly 
behaved  (except  for  the  small  random  perturbations) . 

A  concise  way  to  explain  the  effects  of  the  random  perturbations 
on  a  stable  two  dimensional  manifold  would  be: 

(1)  Any  random  displacement  of  the  system  out  of  the  manifold  will 
vanish,  since  the  system  will  be  drawn  back  into  the  manifold. 

(2)  Any  sequence  of  random  displacements  of  the  system  within  the 
manifold  leads  to  a  two-dimensional  random-walk  motion  within  the 
manifold  superimposed  on  the  motion  due  to  the  vector  field  on  the 


manifold. 
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(3)  Since  the  random  displacements  are  very  small,  the  random-walk 

motion  is  negligible  compared  to  the  motion  induced  by  the  vector 
field. 

It  is  evident,  that  the  above  descriptions  hold  for  the  original  2-torus 
in  the  complete  state  space  as  well  as  for  its  projections  into  the 
individual  three-dimensional  subspaces. 

To  investigate  the  effect  of  random  perturbations  on  the  system 
after  it  underwent  three  bifurcations,  consider  the  motion  of  the  system 
projected  into  one  of  the  three-dimensional  subspaces.  The  motion  will 
be  determined  by  one  of  the  possible  vector  fields  in  the  torus  shell,  as 
described  earlier.  Any  sequence  of  random  perturbations  acting  on  the 
point  will  result  in  a  three-dimensional  random-walk  motion  superimposed 
on  the  vector-field  induced  motion  in  the  torus  shell.  This  can  be 
concluded  by  advancing  the  two-dimensional  case  previously  discussed 
to  three  dimensions.  For  random  displacements  out  of  the  torus  shell 
again  the  considerations  for  the  2-torus  with  the  evident  modifications 
apply. 

One  way  of  looking  at  the  dynamics  within  a  three-dimensional 
subspace  is  to  interpret  the  velocity  components  belonging  to  this  sub¬ 
space  as  time  dependent  variables  of  this  subsystem  under  observation 
and  the  velocity  components  of  the  remaining  system  as  time-dependent 
parameters  of  this  subsystem.  In  other  words,  one  could  construct  a 
non-autonomous  system  of  three  ODE's  for  the  fluid  motion  through  a 
particular  point  in  fluid  space,  if  the  functional  time  dependencies  of 
all  the  other  points  in  fluid  space  are  known.  The  previously  discussed 
structures  of  several  possible  vector  fields  within  the  torus  shell  in 
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the  three-dimensional  space  of  the  three  ODE's  would  be  dependent  on 
the  setting  of  the  other  velocities  in  fluid  space.  That  is,  which  one 
of  the  possible  vector  fields  will  be  realized,  is  determined  by  the 
parameters  of  the  subsystem,  i.e.,  the  velocity  field  in  fluid  space  of 
the  whole  system  minus  the  monitored  subsystem  under  observation.  One 
could  say  therefore,  that  the  structural  stability  of  the  vector  field 
momentarily  realized  in  the  three-dimensional  subspace  depends  on  the 
momentary  setting  of  its  determining  parameters,  i.e.,  the  velocity 
field  in  fluid  space  of  the  whole  system  minus  the  subsystem  under 
observation.  Perturbing  these  parameters  means  perturbing  the  vector 
field  in  the  subsystem  under  observation.  The  main  assumption  of  the 
model  is  now  that  certain  combinations  of  random  perturbations  on  the 
parameters  determining  these  vector  fields  can  change  them  (i.e.,  the 
parameters)  in  a  way  that  the  vector  field  determined  by  them  in  the 
subspace  changes  as  well.  A  sequence  of  these  random  perturbation 
combinations  would  therefore  lead  to  a  sequence  of  random  changes  between 
the  realizable  vector  fields  in  the  subspace,  and  none  of  these  vector 
fields  would  have  a  lifetime  long  enough  to  be  considered  structurally 
stable . 

Mathematical  Formulation 

The  foregoing  considerations  can  also  be  expressed  by  a  mathemat¬ 
ical  formalism.  To  do  so,  the  NS-equation  is  being  written  down  in  its 
most  general  form: 


iXUI  =  X(v(t)) 


(3.1) 


Introducing  a  domain  D  in  (physical)  fluid  space  within  which 
Equation  (3.1)  applies  and  an  infinite-dimensional  (mathematical) 
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state  space  H  for  Equation  (3.1),  V(t)  can  be  interpreted 
physically  as  a  velocity  vector  field  over  D  at  time  t  or 
mathematically  as  a  point  within  H.  This  point  describes  the  state  of 
velocity  of  the  system  within  D  at  time  t.  V(t)  can  therefore  be 
defined  as  the  following  set: 

V(t)  =  (v (r , t) ;  r  e  D},  (3.2) 

with  v(r,t)  as  the  physical  velocity  vector  at  point  r  and  time  t. 

Accordingly,  X  can  be  interpreted  physically  as  the  acceleration 
vector  field  over  D  at  time  t,  or  mathematically  as  a  vector  field 
over  H  describing  the  change  of  the  state  of  velocity  of  the  system 
within  D  at  time  t,  whatever  state  (or  point  in  H)  the  system 
occupies . 

X  as  a  vector  field  over  H  is  time  independent;  all  points  in 
H  can  be  identified  with  all  the  possible  states  of  the  system,  and  a 
change  from  one  particular  system  state  to  another  does  not  depend  on 
a  particular  point  in  time  but  on  the  particular  system  state  only.  A 
particular  point  V  in  H  will  therefore  have  its  time  independent 
"rate-of-change  vector"  X(V),  irrespective  of  the  time  the  state  V 
realizes . 

In  its  physical  interpretation  X  is  time  dependent;  looking  at 
a  particular  point  in  D  relates  to  projecting  X  from  the  infinite¬ 
dimensional  H  into  the  three-dimensional  subspace  H. ,  say,  spanned 
by  the  velocity  components  at  that  particular  point  r^  e  D.  This  loss 
of  dimensions  belonging  to  the  cospace  of  EL  is  automatically  compen¬ 
sated  for  by  introduction  of  time.  One  could  also  say  that  this 
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projection  is  accompanied  with  a  change  of  Eulerian  coordinates  fixed 
within  H  to  Lagrangean  coordinates  travelling  withthe  point 
representing  the  momentary  state  of  the  system.  Everything  time 
independent  within  the  Eulerian  coordinates  has  to  look  time- dependent 
in  the  Lagrangean  coordinates.  X  can  therefore  be  written  as: 

X(V(t))  =  X(t)  =  |3v3(pt:i  ;  r  €  d]-  (3.3) 

(physical;  Lagrangean  coordinates) 

or 

X(v(t))  ;  V(t)  £  h},  (3.4) 

(mathematical;  Eulerian  coordinates) 

whereby  the  second  formulation  (3.4)  is  the  more  general  one 
relating  to  Equation  (3.1). 

Geometrically,  the  solution  space  H  can  be  understood  as  the 
cartesian  product  of  the  three-dimensional  orthogonal  subspaces  EL  ; 
one  for  each  point  r^  e  D  . 

★ 

H  =  x  H.  ,  H.  =  R3  ,  v(r.,t)  e  H.  (3.5) 

„  l  l  l  l 

r .  eD 

l 

This  splitting  carries  over  to  the  vector  field  X  as  a  direct 
sum  over  all  r.  e  D  : 

l 

*  3 

This  should  be  read  as:  "H.  is  a  copy  of  ]R  i.e.,  transitivity 

3  1 *  3 

does  not  apply:  EL  =  ER  ,  EL  =  1R  ,  but:  EL  ±  EL  ,  for  i  f  j. 
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X(V(t))  =  ©  X  (V(t))  (3.6) 

r.eD 

1 

with  X^(V(t))  as  a  three-dimensional  vector  at  point  V(t)  in  H 
(not  in  H^!!)  with  its  components  only  in  the  subspace  hL  ,  or,  in 
the  physical  picture,  the  acceleration  vector  at  point  r^  ,  "velocity 
state"  V,  and  time  t.  Equations  (3.2),  (3.3)  and  (3.6)  show  now 
that  Equation  (3.1)  can  be  decomposed  accordingly  into  terms  like: 

9v(r. ,t) 

- ^ -  =  X.(V(t))  (3.7) 

In  other  words.  Equation  (3.1)  can  be  represented  as  a  system 
of  equations  of  type  (3.7);  one  for  each  r^  £  D.  The  further  formal 
"physical  decomposition"  of  Equation  (3.7)  into  three  coupled  scalar 
equations  for  the  velocity  and  acceleration  components  at  r^  is  now 
evident: 

3vk 

TT  =  ;  k  =  *  (3*8) 

Each  subsystem  of  the  type  (3.7)  now  depends  not  only  on  the 
velocity  v(r.,t)  of  the  point  r^  ,  to  which  it  applies,  but  (as  can 
be  seen  by  the  definition  in  Equation  (3.2))  on  all  the  velocities  of 
the  other  points  r  £  D  as  well,  i.e.: 

3v(r.  ,t)  ^ 

- ^ -  =  X.(v(r.,t);  V (t)  \i)  (3.9) 


with: 
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V(t)  \  i  =  |  v(r,t) ;  r  'e  D  -  {r^ 


(3.10) 


The  right  hand  side  of  (3.9)  can  now  be  written  as: 


3v(ri,t) 

at 


=  x. 


V(t)\i(v(ri>t)) 


(3.11) 


to  indicate  the  role  of  the  elements  of  V(t)  \  i  as  (time  dependent) 
parameters.  ^y(t)\i  a  three-dimensional  vector  field  in  FL  . 
Physically,  however,  it  is  identical  to  X^(V(t)),  since  it  performs 
the  same  physical  function. 

To  investigate  the  behaviour  of  the  described  system  in  the 
presence  of  a  torus  as  stable  invariant  manifold,  some  geometric 
properties  of  tori  embedded  in  higher  dimensional  spaces  are  being 
stated.  Looking  at  a  circle  embedded  in  such  a  space  in  a  way  that 
its  plane  or  symmetry  axis  is  not  parallel  to  any  of  the  coordinates 
spanning  the  space,  it  is  evident  that  if  two  points  on  this  circle 
have  more  than  one  identical  coordinate,  the  two  points  have  to  be 
identical.  Moreover,  it  is  evident  that  amongst  the  set  of  possible 
cartesian  coordinate  systems  for  this  higher  dimensional  space,  the 
subset  with  coordinates  parallel  to  the  circle  plane  or  symmetry  axis 
has  Lebesgue  measure  zero.  For  a  2-torus,  the  argument  holds  as  well, 
except  that  more  than  two  coordinates  have  to  be  identical  for  the  rest 
of  the  coordinates  to  be  identical  as  well.  This  can  be  stated  for 
higher-dimensional  tori  as  well  in  the  following  general  way: 

Let  M  be  a  manifold  embedded  in  ]Rn  and  Pn  ,  P9  e  M. 

J.  ^ 


Pn  =  (x  }n 
1  ll  i=l  * 


n 


P9  =  {x.9}.  , 

2  i2  i=l 


(3.12) 
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since  they  are  points  in  ]Rn  as  well. 

If  there  exist  subsequences  S™  ,  S™ 

of  Pn  Pn  * 
or  r1  ,  r2  . 


within  the  coordinate  strings 


cm  j-  -)  m 

si  - 

S™  =  (x. 
2  J 

m 

2  j  =  1 

9  m  < 

n,  -  x.2Vj 

(3.13) 

if: 

m  >  1  and 

M  =  T1 

then 

Xil  = 

x12Vi  j 

m  >  2  and 

2 

M  =  T 

then 

Xil  = 

Xi2Vi 

(3.14) 

m  >  k  and 

V 

M  =  T 

then 

Xil  = 

xi2Vi  ) 

except  for  a  set  of  cartesian  coordinate  systems  of  Lebesgue  measure 
zero. 

Assuming  the  existence  of  a  circle  or  a  2-torus  as  invariant 
manifold  M  in  H  and  two  points  V(t^),  V(t9)  on  M: 

V(t  ),  V(t2)  e  M  c  H;  M  =  T1  or  T2  (3.15) 


such  that  for  the  restriction  or  projection  IL  of  V(t^),  V(t©  into 
one  of  the  three-dimensional  subspaces  PL  equality  holds,  i.e.: 

ni(V(t1))  =  vCr.,^)  =  n.(V(t2))  =  v(r.,t2)  (3.16) 

(I)  (II)  (HI)  (IV) 


and: 


X(M)  =  ©  X. (M)  (Whitney  Sum) 

r.eD  1 

l 


(3.16a) 
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one  can  conclude  by  the  foregoing  argument.  Equation  (3.16)  (II)  =  (IV) 
and  by  identifying  the  three  components  of  v(r^,t^),  v(r  ,t^)  with 
the  members  of  the  subsequences  S™  ,  S™  that: 

V(tx)  =  V(t2)  (3.17) 


or: 


VpXi^i’V5  '  ?V(t2)\i(7(ri't2))’  M  =  ?1  °r  ?2  (3-18) 


That  means,  one  and  the  same  point  v(r^)  fixed  within  FL  holds 

— 14  —S. 

one  and  the  same  field  vector  X,.  -  v  .  (v(r . ) )  for  all  times.  In 

V(t)\i  v  l 

other  words,  if  the  invariant  manifold  M  c  H  is  a  circle  or  2-torus, 
the  vector  field  X,..  . *  .(v(r.,t))  associated  with  it  in  H.  will  be 
independent  of  V(t)  \i,  V(t)  and  thereby  independent  of  time  as  well, 
If  the  invariant  manifold  M  c  H  is  a  3-torus  one  cannot  conclude 
Equation  (3.17).  Since  only  three  coordinates  have  to  be  identical 
to  satisfy  Equation  (3.16),  m  =  k  =  3  holds.  Therefore,  condition 
(3.14)  is  not  fulfilled,  or: 


XV(t1)\ifV^ri,tl^  ^  XV(t2)\itV(;ri,t2:)')  *  M  T 


(3.19) 


is  possible,  even  if: 

v(ri,t1)  =  v(r.,t2),  tx  ^  t2  (3.19a) 

Therefore,  the  vector  field  in  H.  associated  with  an  invariant 
’  i 

3 

manifold  M  =  T  c  H  can  change  with  V(t)\i  =  (v(r,t);  reD  -  {r.}}, 
that  is,  with  all  the  velocities  at  all  points  r  e  D,  except  r^  . 
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One  can  now  introduce  the  small  random  perturbation  distribution 
VR(t)  over  D  at  time  t  with  a  negligibly  small  rise  time:* 

VR(t)  =  (vR(r,t);  r  e  D}  (3.20) 

and  impose  its  reduced  form: 

VR(t)  \  i  =  {vR(r,t);  r  e  D  -  {r^j-  (3.21) 


upon  V  (t)  \  i  : 


V(t)  \  i  +  VR(t)  \  i  =  { v(r,t)  +  vR(r,t)  ;  r  e  D  -  { r±}  f  (3.22) 


This  leads  to: 


XV(t)\iCv(ri’t):>  +  XV(t)\i+VD(t)\i(V(ri't)) 

K 


(3.23) 


Assume  now  for  certain  perturbation  distributions  VD  (t  )  the  following: 

K  1 

V(tp  \i  +  VR(tp  \  i  =  V(t2)  \  i,  (3.24) 

and : 

v(ri,t1)  =  v(ri,t2),  t1  t  t2  (3.24a) 

1  2 

Equation  (3.18)  show  now  that  for  M  =  T  or  T  the  vector  field  in 
EL  remains  the  same  or  stable  under  VR(t)  or  VR(t)  \  i  respectively, 

3 

but  for  M  =  T  inequality  (3.19)  shows  that  the  vector  field  can 
adopt  a  structure  assigned  to  a  completely  different  point  in  time,  i.e., 
t ^ ,  since  condition  (3.24a)  is  the  same  as  condition  (3.19a),  which 

3 

permits  in  Equation  (3.19)  for  M  =  T  .  Therefore,  the  following 
transition  occurs: 

•k 

The  time  it  takes  for  the  perturbation  to  reach  its  maximum. 
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Vpxi^i’V5  *  V2)siCv(ri.t2»  -  ?vCt JxiW'i'1!’  (3-25^ 


The  equation  holds  since  vfr^t^)  =  v(r^,tp,  (more  realistic: 
v(r. ,t„)  C5iv(r.,t,),  since  the  rise  time  for  V„(t)  is  finite 
although  negligibly  small) .  This  means  that,  although  the  individual 
random  perturbation  on  the  individual  point  in  D  might  be 
infinitesimally  small,  certain  distributions  of  them  (compatible  with 
Equation  (3.24))  imposed  on  a  set  of  points  in  D  will  uncouple  a  point 
r^  e  D  for  the  rise  time  of  the  perturbation  distribution  from  the  rest 
of  the  vector  field  over  D  by  breaking  down  the  vector  field 
X  )\i(v(ri,t))  i-n  Hq  and  then  recouple  it  by  realizing  a  new 
field  -j y^(v(r^,t))  almost  randomly  selected  from  a  set  {Xy^^}, 

say,  permissible  by  Equation  (3.11).  This  new  vector  field 
Xy(t  -j\^(v(r^,t))  might  not  necessarily  be  C  equivalent  to 

X..,  .  v  .  (v  (r .  ,  t) )  ,  which  would  mean  that  some  of  the  elements  X,r,  .  >  . 

V(t^) \i  v  i’  V (t) \i 

of  the  set  ^y(t)\i^  are  not  structuraHy  stable.  Clearly,  a 
sequence  of  such  random  distributions  leads  to  chaotic  behaviour  at  the 
individual  point  r^  e  D  and  this  therefore  generalized  directly  to 
the  entire  D. 

To  complete  this  description,  some  consideration  is  given  to  the 
already  mentioned  set  of  cartesian  coordinate  systems  of  Lebesgue 
measure  zero,  for  which  the  relations  (3.14)  do  not  hold.  That  is, 
if  the  invariant  manifold  is  a  k-torus,  two  different  points  on  it  can 
be  in  more  than  k  coordinates  identical.  The  random  perturbed 
behaviour  of  the  system  on  a  3-torus  will  essentially  be  the  same  as 
already  described,  except  for  a  possible  change  in  the  set  of 
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selectable  vector  fields  (X...  .  .}  on  H.  .  For  a  circle  or  2-torus, 

V(t)\i  1 

however,  for  some  set  L,  say,  of  subspaces  FL  (with  Lebesgue  measure 

zero),  inequality  (3.19)  instead  of  Equation  (3.18)  will  hold.  That 
means,  for  one  or  two  bifurcations  some  points  r^  (associated  with 
the  FL  €  L)  can  show  chaotic  behaviour.  At  first  glance,  it  appears 
as  if  this  kind  of  behaviour  can  be  induced  deliberately  by  the  way  the 
coordinate  system  spanning  H  is  being  set  up  or  rotated.  One  might 
be  able  to  select  one  belonging  to  the  set  (of  Lebesgue  measure  zero) 
for  which  the  relations  (3.14)  do  not  hold  at  will.  However,  H  is 
spanned  by  the  (see  Equation  (3.5))  and  the  structure  of  any 

combination  of  them  is  determined  by  the  structure  of  the  r^  e  D , 
i.e.,  the  geometry  of  D.  The  three-dimensional  cartesian  coordinate 
systems  spanning  the  individual  ' s  can  be  rotated  in  any  way  with¬ 
out  changing  the  applicability  or  inapplicability  of  the  relations 
(3.14)  for  the  individual  r^'s.  To  obtain  such  an  effect,  the 
individual  H^'s  would  have  to  be  spanned  by  velocity  components 
belonging  to  different  points  r  6  D  and  could  therefore  not  be 
associated  with  one  particular  r^,  but  two  or  three.  Obviously,  this 
would  be  meaningless  within  this  model. 
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PART  2 


Technical  Realization  and  Approximation  of  the  Model 

As  already  mentioned  in  the  introduction,  there  are  two  levels 
of  analysis  of  the  fluid  dynamic  equations.  One  level,  the  qualitative 
analysis  and  description  of  the  mechanism  for  turbulence  was  the  subject 
of  the  previous  chapter.  In  the  present  part,  the  approach  is  more 
application-oriented  in  that  a  technical  basis  for  computer  simulation 
is  developed.  Before  the  main  problem  of  discretization  of  the  fluid 
domain  D  is  addressed,  a  particular  type  of  NS-equation  is  derived, 
which  permits  one  to  keep  a  certain  generality  w.r.t.  the  selection 
of  various  flow  structures. 

4.  THE  NS-EQUATIONS  FOR  PERTURBATIONS 

For  a  velocity  vector  U  at  point  r  and  time  t,  the  NS- 
equation  is  written  as: 

V 

9U(rt)  =  _u(r>t)[v  ;  U(r,t)]  +  vAU(r,t)  -  Vp(r,t)  (4.1) 

O  L 


with 


as  the  advective  term,  and 

v  =  kinematic  viscosity, 
p(r,t)  =  pressure  at  point  r  and  time  t. 


(4.2) 


35 


36 


For  convenience,  an  advective  operator  A  can  be  introduced: 


A[ab]  =  a[V  :  b]  +  b[V  :  a] 


(4.3) 


Clearly,  this  operator  is  symmetric  and  bilinear  in  a  and  b. 
Equation  (4.1)  becomes  now 


9U(r,t)  _ 
at 


=  -  4  [U (r , t) U(r , t)  ]  +  vAU(r,t)  -  vp(r,t) 


(4.4) 


and 


holds  also,  if  U  develops  a  perturbation  : 


U(r,t)  +  U(r,t)  +  Vx(r,t) 


3  (U+Vl)  A 


—  =  -  j  [(u+vpai+vp]  +  vA(U+V1)  -  vp 


(4.5) 


(4.6) 


(For  reasons  of  conciseness,  the  dependencies  (r,t)  are  not 
written  out  in  Equation  (4.6)). 

Since  A  is  bilinear.  Equation  (4.4)  can  be  subtracted  from 
Equation  (4.6)  to  obtain  an  equation  for  the  perturbation  : 


3  V 


±  =  -A[U  vp  -  |  [VjVj]  *  vAVj 


(4.7) 


_ -X  _ 3k 

Using  the  same  method  for  a  perturbation  V0  on  top  of 


V1  ^  V1  +  V2 


(4.8) 


leads  to  an  equation  for  V0  : 


3y2 

IT  -  -A[UV2]-  AtVjVj] 


2  [V2V2]  +  vAV2 


(4.9) 
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A  third  perturbation  V 3 


on  top  of  ^2  • 


-> 


(4.10) 


leads  to  an  equation  for 


=  -A[UV3]  -  A[V1V3]  -  A[V2V3]  -  J  [V3V3]  +  vAV: 


(4.11) 


Naturally,  these  equations  for  the  perturbations  can  also  be  obtained 
without  the  introduction  of  the  advective  operator  A.  Their 
derivation,  however  would  be  rather  lengthy. 

A  more  accurate  and  realistic  approach  would  require  the 
introduction  of  a  pressure  perturbation  as  well,  e.g.: 

P  P  +  V1  (4.12) 


However,  this  would  require  an  evolution  equation  for  the  pres¬ 
sure.  Unfortunately,  such  an  equation  does  not  exist.  Moreover,  it  is 
evident  that  this  neglect  is  legitimate  for  the  assumption  of 
isotropic  pressure  perturbations,  since  their  gradient  vanishes. 

Equations  (4.7),  (4.9),  (4.11)  determine  the  stability  of  a 
given  solution  by  solving  for  the  time  development  of  a  perturbation. 
Assume  a  given  flow  field  U(r,t)  as  a  solution  of  Equation  (4.4)  and 
substitute  it  into  Equation  (4.7).  The  perturbation  V  (r,t)  will 
vanish  for  stable  U(r,t)  and  grow  unbounded  or  become  periodic  for 
U(r,t)  unstable  or  undergoing  a  Hopf  bifurcation.  U(r,t)  +  V  (r,t) 
can  now  be  substituted  into  Equation  (4.9)  to  obtain  the  perturbation 
V2(r,t)  and  U(r,t)  +  V^(r,t)  +  V2(r,t)  in  Equation  (4.11)  works 
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accordingly  with  V^(r,t). 

This  particular  approach  was  motivated  by  the  availability  of 
a  FORTRAN  program  which  investigates  (Hopf)  bifurcational  behaviour 
of  ODE's.  The  program  can  be  used  on  a  discretized  version  of 
Equation  (4.7),  i.e.,  a  system  of  ODE's,  to  test  a  given  stationary 
(discretized)  flow  field  U(r)  for  Hopf  bifurcations.  If  a  Hopf 
bifurcation  occurs,  U(r)  and  the  periodic  solution  V  (r,t)  can  be 
tested  in  Equation  (4.9)  (discretized)  for  a  second  Hopf  bifurcation 
which  leads  to  a  2-torus  as  invariant  manifold  in  solution  space. 

After  that.  Equation  (4.11)  can  be  used  to  test  for  the  third  Hopf 
bifurcation.  Evidently,  this  method  "lives"  on  the  particular 
structure  of  the  advective  terms  in  the  NS-equation  or  Equations  (4.7), 
(4.9),  (4.11)  respectively.  It  is  this  term  which  permits  the 
construction  of  evolution  equations  for  perturbations  without  elimi¬ 
nating  the  underlying  flow  out  of  the  equation. 

However,  the  nonautonomous  nature  of  Equations  (4.9)  and  (4.11) 
causes  problems  for  the  application  of  the  bifurcation  program,  since 
time-dependent  functions  instead  of  stationary  distributions  have  to 
be  tested  for  bifurcations.*  An  attempt  to  deal  with  this  problem  is 
described  later  in  the  chapter. 

Discretization  of  the  NS-equations** 

To  obtain  numerical  approximate  solutions  for  equations  of  the 
type  (4.7),  (4.9),  (4.11),  they  have  to  be  approximated  by  a  system 

■ k 

See  also  Chapter  2 

• k  ★ 

See  Appendix  D  for  an  introduction  to  previous  works. 
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of  ODE’s.  Such  a  system  is  obtained  by  selecting  a  discrete  finite 
set  S,  say,  of  points  rs  from  the  fluid  domain  D  and  constructing 
a  subsystem  of  three  ODE's  for  the  three  components  of  the  unknown 
perturbation  velocities  V(rs)  at  each  rg  e  S.  These  different 
subsystems  have  to  be  coupled  mutually  to  account  for  the  interaction 
between  the  different  r^  e  S.  In  the  next  section  such  a  system  of 
ODE’s  for  Equation  (4.7)  is  derived  using  finite  element  techniques 
and  the  Galerkin  method.  At  the  same  time,  some  of  the  theory  and 
motivation  behind  this  approach  is  explained.  For  instance,  it  is 
being  shown  how  incompressibility  and  Dirichlet  boundary  conditions 
can  be  obtained  by  appropriate  construction  of  basis  functions. 
Furthermore  an  expansion  from  the  ODE-system  for  Equation  (4.7)  to 
systems  for  Equations  (4.9)  and  (4.11)  is  performed. 

5.  THE  GALERKIN  METHOD 

One  of  the  most  widely  used  and  successful  discretization- 
approximation  techniques  is  the  Galerkin  method.  It  is  applied  here 
in  connection  with  some  kind  of  function  space  representation,  i.e., 
the  functions  U,  appearing  in  the  Equation  (4.7)  are  being 
represented  as  vectors  in  the  infinite- dimensional  function  space  H 
introduced  in  Chapter  3: 

A  |vx>=  -  I  U[V  :  Vx]  >  -  I  Vx[v  :  U]  >  +  v|  AV:  >  -  |  [V  :  V]  >*  (5.1) 

■ k 

For  reasons  of  conciseness  the  spatial  and  temporal  dependencies 
(r,t)  are  not  written  out. 


"C 


40 


This  is  Equation  (4.7)  written  in  function-space  notation  without  the 
use  of  the  advective  operator  A.  The  symbol  |  >  is  borrowed  from 
the  formalism  used  in  quantum  theory  and  denotes  a  point  or  vector  in 
function  space,  i.e.,  in  H  here. 

The  '’subdecomposition"  into  fluid  space  components  is  also 
evident,  e.g.  : 


(5.2) 


Again  following  the  quantum  mechanical  formalism,  the  identity 
operator  I  is  introduced: 


I 


l 

reD 


A  >  <  A 

r  r 


(5.3) 


A^  is  the  three-dimensional  Dirac  delta  function  of  the  point 

r  e  D  (see  Appendix  B) .  The  evident  interpretation  of  definition 

(5.3)  within  the  concept  of  function  spaces  or  infinite-dimensional 

vector  spaces  is  that  of  a  point  or  vector  |  v>  ,  say,  being 

decomposed  by  I  into  the  components  parallel  to  the  coordinates  of 

the  infinite-dimensional  system  defined  by  the  orthogonal  vectors 

I A  >  ,  i.e., 
r 


1 1  v  >  =  i 

reD 


A  >  <  A  V> 

r  r  1 


(5.4) 


The  <A^  I  V>  can  be  understood  as  the  coordinates  or  component  values 
of  |  V>  within  the  system  spanned  by  the  |  A^  >  ,  r  e  D,  i.e.,  the 
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scalar  products  between  the  |  >  and  |  V>  .  * 

In  order  to  separate  the  perturbation  velocities  V  out  to 

the  right  the  identity  operator  I  will  be  applied  to  the  left  of 

each  of  these  perturbations  in  Equation  (5.1).  After  that,  a  finite 

set  S  of  points  r^  is  selected  from  the  entire  domain  D.  This 

restricts  the  calculations  of  the  perturbation  field  V^(r,t),  r  e  D 

to  the  calculations  of  the  discretely  distributed  V, (r  ,t) ,  r  e  S: 

Inn 

and  reduces  the  partial  differential  NS-equations  (4.7)  to  a  system 
3  x  n  ODE's,  where  N  is  the  cardinality  of  S.  (This  includes  the 
fluid  space  subdecomposition  according  to  Equation  (5.2)  as  well.) 

The  spatial  derivatives  introduced  therein  by  the  Del  and  Laplace 
operators  could  now  be  replaced  by  finite -difference  approximations 
constructed  of  elements  like  e.g.: 


9f (r  )  f(r  . ) -f (r  J 
n  n+1  n-1 

3r  ~  2h 


(5.5) 


with : 


r^  =  selected  point 

r^+^  =  selected  points  in  neighbourhood  of  r^ 

h  =  spatial  distance  between  r  ,  r  .  or  r  ,  r  , 
r  n5  n+1  n  n-1 


resp . ** 


This  and  various  other  finite- difference  approximations  will  not  pre- 

_S» 

serve  the  incompressibility  condition  div  U  =  0  of  an  initial 
velocity  field  U  in  the  subsequent  numerical  iterations,  regardless  of 


Introduction  of  the  duality  concept  between  bra  (<|)  and  ket  (|  >) 
vectors  is  not  necessary,  since  only  real-valued  functions  are  used 
in  this  context. 

Naturally  a  discretization  with  equispacing  is  assumed. 
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whether  they  develop  the  evolution  NS-equation  in  time  or  just  iterate 
a  steady  state  NS-equation  via  algorithms  towards  a  higher  accuracy. 

One  therefore  depends  on  a  method  that  preserves  incompressibility  . 

Such  a  method  is  most  efficiently  dealt  with  and  explained  in  the 
function  or  state  space  context. 

The  state  space  of  the  NS-equations  is  a  subspace  of  the 
function  space  represented  by  the  Dirac  functions  as  basis  functions 
,  r  e  D.  That  is,  application  of  the  identity  operator  I  defined 
in  Equation  (5.3)  will  represent  a  NS-solution  as  a  point  in  a  space 
spanned  by  the  vectors  |a^>,  r  e  D.  Within  this  state  space  there 
is  a  subspace  V  in  which  the  condition  div  U  =  0,  V|U>e  V  holds. 
The  obvious  choice  for  a  set  of  basis  functions  spanning  such  a  space 
is  a  set  of  scalar  functions  W  (y),  x,  y  e  D  which,  when  multiplied 

X 

with  a  velocity  field  U  over  D  yield  a  divergence-free  vector  field, 
i .  e.  : 

div(Wx(y)  •  U(x))  =  0,  (5.6) 

where  the  derivatives  are  taken  w.r.t.  y  and  the  velocity  U(x)  is 
the  value  of  the  field  taken  at  the  point  x  ;  that  is,  x  is  not  to 
be  understood  as  a  variable  in  Equation  (5.6).  (It  is  evident  that 
these  basis  functions  don't  have  to  be  orthogonal  to  each  other  such 
as  the  Dirac  delta  functions  .  Their  supports  within  D  may 
overlap . ) 

Replacing  each  Dirac-type  basis  function  A  ,  x  e  D  by  W  (y) , 
x,  y  e  D  leads  to  the  following  representation  of  the  solution  U  : 


43 


l  I  W  (y)XW  (y)|  U>  =  l  W  (y)  /  W(y)U(y)dy 
xeD  X  X  x«D  D  X 

=  l  W(y}U,  (5.7} 

xeD 

i.e.,  a  linear  combination  of  an  infinite  number  of  basis  functions 

v  _ ^ 

W  (y)  with  some  type  of  mean  values  U  of  U  as  coefficients. 

If  the  solution  U  had  nonzero  divergence,  application 
of  Equation  (5.7)  would  project  U  into  a  divergence-free  or 
solenoidal  subspace  of  the  solution  space.  Without  knowledge  of  the 
entire  velocity  field  U  ,  the  means  U  cannot  be  calculated  and 
therefore,  only  the  point  values  U(r  )  at  the  points  r^  e  S 
selected  from  D  are  being  used.  This  discretization  produces  the 
following  approximation  to  Equation  (5.7): 


r 

n 


I 

eS 


(y)><  Ar 

n 


U> 


r 

n 


l  (y)  /  Ar  U(y)dy 

eS  n  D  n 


=  I 

Vs 


(y)U(rn) 


(5.8) 


Application  of  the  herein  defined  "discretization-approximation- 
projection"  operator 


dap 


n 


I 

eS 


t*  (y)>  <  A 


n 


n 


(5.9) 


to  the  perturbation  velocities  in  Equation  (5.1)  leads  to  a 

system  where  the  spatial  derivatives  are  taken  analytically  on  the 
's  instead  of  using  some  finite  difference  approximation  over 


n 

discrete  points  as  shown  in  Equation  (5.5) 


Applying  therefore  on 


IV 


in  Equation  (5.1)  from  the 
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left  gives: 


at 


l  |wr  (y)XAr  |v  >  =  -  l  |iT[V:Wr  (y)]XAr  |vx> 

r  eS  n  n  r  eS  n  n 

n 


r  eS 
n 


l  I  [V  :  U]Wr  (y)X4r  |v  > 

r  eS  n  n 

n 


+  v  l  I A  •  W  (y)X  A  |v,  > 

r  eS  rn  n  1 

n 


I  ]wr  (y)VWr  (y)XA  |vi><Ar  |V  >  (5.10) 

r  ,r  eS  n  m  n  m 

n  m 


here: 


Clearly,  this  follows  directly  from  Equation  (4.2).  The 
nonlinear  term  is  derived  as  follows: 


(5.10a) 


(5.10b) 
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Vx  [V  :  V1]> 


'dap’ 


l  |wr  (y)><q.  |Vj>  l  lvwr  (y)Xar  |vx> 

r  eS  n  n  r  eS  m  m 

n 


r  eS  m 
m 


r  ,r  eS  n 
n  m 


(y)><A  V  >  vw  (y)><A  .> 

r  1  r  r  1  1 

n  mm 


I  |wr  (y)vwr  (y)XAr  | v  > < ar  |v  > 

r  ,r  eS  n  m  n  m 

n  m 


Next,  Equation  (5.10)  is  multiplied  from  the  left  with  the  operators: 


<wr  (y) 

o 


1  0  0 

0  10  =<W  (y)  1 1? > 

0  0  1  o 


r  €  S 
o 


(5.11) 


to  obtain  a  system  of  3  x  N  ODE's: 


3t  q<Wr  W><Ar  I  V1  > 

r  eS  o  n  n 

n 

=  -  l  <K  Cy)I3|U[V:Wr  (y)]Xar  |v1> 


r  eS  o 
n 


n 


n 


l  <>*r  (y)i3l  [v  :  u]wr  ty)Xar  |v  > 

r  sS  o  n  n 

n 


+  V  l  <Wr  Cy)  Ij  |  AWp  (y)XA  |v  > 

r  eS  o  n  n 

n 


l  Cy)I3|Wr  (y)VWr  (y)XAr  [V1><Ar Jv^  ;  roe  S. 

(5.12) 


r  ,r  eS  o 
n  m 


n  m 


n 


m 


This  is  a  system  of  known  coefficient  matrices  and  vectors 
consisting  of  the  unknown  perturbation  velocities  at  the  points 

r^  £  S.  The  3x3  unit  matrix  I  controls  proper  subdecomposition 


* 

The  arrow  translates  as: 


"application  of  yields". 
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of  the  coefficients  into  their  fluid  space  components. 

To  show  the  structure  of  these  matrix  elements  or  coefficients, 
especially  in  their  fluid-space  subdecomposition,  they  are  being 
written  out  explicitly  as  integrals  and  at  the  same  time  technical 
abbreviations  are  introduced  below. 


VE  <Wr  Wjl"!  (y)>  =  /  Wr  Cy)Wr  (y)dy 
o  n  Don 


1  0  0 
0  1  0 
0  0  1 


(5.13) 


WUV  =  <W  (y)I_|U[v"  :  W  (y)  ]  > 
on  r  7  31  L  r  7  J 

o  n 


/  wr  (y) 

D  o 


ux(y) 


8Wr  (y) 
n 


3Wr  (y)  3Wr  (y) 

+  My)  TT7 - +  U,(y)  - 


y2 


y3 


'l 

0 

o' 

dy 

0 

1 

0 

.0 

0 

1 . 

(5.14) 


,  V:  U 


W  E  <  W  (y  )  I  _  [V:  U]W  (y)> 
on  r  7  3)  J  r  7 

o  n 


pu  1 

"bu  I 

/  wr 

D  o 

Cy)w 

n 

(y) 

[*;] 

dy 

f  W 

D  ro 

Cy)wr 

n 

(y) 

k 

dy 

/  W 

7  r 

D  o 

(yjwr 

n 

(y) 

l 

_3y  3_ 

dy 

[bu  1 

rsu  - 

"bu  " 

f  w 

D  ro 

Cy)wr 

n 

(y) 

dy 

/  \ 

D  0 

Cy)w 

n 

(y) 

z 

_3y2  j 

dy 

f  W 
j  r 

D  o 

(y)wr 

n 

(y) 

Z 

_3y  3_ 

dy 

"BU/ 

rsu  n 

rau  i 

f  w 

D  ro 

(y)w 

n 

(y) 

5 

-8yl. 

dy 

f  w 

J  r 

D  o 

(y)Wr 

n 

(y) 

3 

_3y2_ 

dy 

/  W 

7  r 

D  o 

(y)Wp 

n 

(y) 

5 

>3. 

dy 

/  W  (y)W  (y) 
D  o  n 


BU, 


By 


1 


BUX  3U 


By:  9y  2  3y: 

^2 

3yi 
3U, 


3U  8U2 

3y9  3y? 


^3  3_h 

3y2  8y3 


dy 


(5.15) 
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W 


A  _ 


on 


E  < W  (y)lJAW  (y )  > 


o 


n 


=  /  w  (y) 


D  ro 


a2wr  (y)  a2wr  (y)  a2wr  (y) 

n  n  n 

- ^ -  +  - - -  +  - - - 


8y 


1 


3y: 


3y- 


"1 

0 

o" 

dy 

0 

1 

0 

_0 

0 

1_ 

(5.16) 


For  the  matrix  elements  relating  to  the  nonlinear  term 


W  e  <W 
onm  r 


Cy)  1-7  I  W  (y)  vw  (y)> 


n 


m 


(5.17) 


a  more  detailed  analysis  is  appropriate.  On  function  space  level  the 
term  is  written  as: 


l  <Wr  (y)I3|W  (y)VWr  (y)X4r  1^X4  |  Vj  > 

r,rgSo  n  m-  n  m 

n  m 


The  terms  <W_^  (y)  I  |  W  (y)VW_^  (y)  >  are  elements  of  a  N  x  N‘ 
o  n  m 

(o  counts  rows  and  n  •  m  counts  columns) ,  whereas  the 


matrix 


<Ar  | Vx> <  A  | >  are  components  of  a  N  -dimensional  column  vector, 
n  m 


Since  <A  IV.XA  I V,  >  =  <  A  IV.XA  I V.  >  these  N  columns 


1 


1 


1 


1 


n  m  m  n 

and  vector  components  can  be  reduced  to  N(N+l)/2  columns  or  vector 
components  respectively. 

On  fluid- space  level  one  has: 
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Since  each  coefficient  matrix  element  of  the  linear  terms  decom¬ 
poses  into  a  3x3  matrix  on  fluid-space  level,  one  actually  has  a 
3N  x  3N  coefficient  matrix  for  each  linear  term.  To  multiply  the 
inverse  of  the  3N  x  3N  coefficient  matrix  of  the  left-hand  side  of 
Equation  (5.12)  with  the  right  hand  side,  one  needs  for  the  nonlinear 
terms  just  as  for  the  linear  terms  a  coefficient  matrix  with  3N  rows. 
Therefore  the  three-dimensional  coefficient  row  vector  of  the  above 
expression  has  to  be  transformed  into  a  matrix  with  three  rows,  and  the 
3x3  matrix  of  the  products  of  the  perturbation  velocity  components 
into  a  column  vector  of  the  same  dimension  as  columns  in  the  matrix- 
transformed  3-dimensional  row  coefficient  vector.  The  definition  of  a 
dyadic  product  between  the  3x3  unit  matrix  I  and  the  3-dimensional 
row  coefficient  vector  casts  the  situation  of  fluid  space  level  into  a 
3x9  coefficient  matrix  and  a  9-dimensional  column  vector  consisting 
of  the  3  columns  of  the  3x3  product  matrix: 


1  0  0 
0  10; 
0  0  1 


3Wr  (y)  3Wr  (y) 

Wr  Cy)Wr  Cy)  “  3y.  ’  Wr  (y)Wr  W  '  37; 

on  1  o  n  2 


Cy) 


Wr  Cy)  — 

o  n 


m 


y- 


V, ,  (r  )V  (r  )  V  (r  )V10(r  )  Vi  (r  )V1_(r  ) 

11  llv  m  11  ^  n  12K  nr  11^  nJ  13^  nr 

V  (r  )V  .  (r  )  V  _(r  )V  0(r  )  V  _(r  )V  _  (r  ) 

12  n  11  m  12v  n  12  m  12  n  13  m 

,  V._(r  jV^fr  )  V  (r  )V  (r  )  V,_(r  )V,„(r  ) 
L  13v  n  11  m  13  n  12v  m  13  n  lo  m 


v 
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ii 


(5.12a) 
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The  nonlinear  coefficient  matrix  (cm)  elements  appear  now  as: 


onm 


<Wr  Cy)I3|Wr  (y)VWr  (y) > 
o  n  m 


=  /  w  (y)W  (y)(3  X  9  cm)  dy 
D  o  n 


(5.13) 


With  the  introduced  abbreviations  the  system  (5.12)  may  be  written  as: 


y  W  4-  V,  (r  ) 

L0  on  9t  lv  n^ 
r  eS 
n 


r  /  nv->.  v  •  u-*-  r  v  \ 

=  y  -Wu  Vn  (r  )  -  W  ^  (r  )  +  vW  V,  (r  )  -  Y  W  V.  (r  )V.  (r  ) 

L  0  \  on  1  n  on  I  n  on  1  n  onm  I  n  1  m  l 

r  eS  \  r  eS  / 

n  m 


r  eS 
n 


,  7UV  ,  rV :  u  ,  ,A 
-W  -  W  +  vW 
on  on  on 


I  WV  V,  (r  ) 
u  onm  1  m 


r  eS 
m 


V.  (r  )  ;  r  e  S 
I  n  o 


(5.14) 


As  can  be  easily  seen  by  comparing  the  nonlinear  term  in  Equation  (5.14) 
with  the  related  term  in  Equation  (4.7)  (using  the  definition  of  the 
advective  operator  (4.3)),  the  discretization  performs  the  following 
transition: 


vpv  :  vp 


discretization 


r  ,r  eS 
n  m 


W  V,  (r  ) V  (r  ) ;  r  e  S  (5.15) 
onm  1 ^  n  1  mJ  o 


The  discretization  procedure  shows  that  this  transition  can  be 
generalized  directly  to: 


V.  [V  :  V.] 
i  3 


discretization 


r  ,r  eS 
n  m 


W  V.  (r  )  V .  (r  )  ;  r 
onm  i  n  j  m  o 


e  S  (5.16) 


or: 


v- 
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V • [V  :  V.  ] 
3  i 


discretization 

- y 


(5.17) 


where  the  indices  n,m  in  W  change  places.  This  difference 
between  (5.16)  and  (5.17)  follows  easily  from  the  definition 
(5.13)  . 

With  the  transitions  (5.16),  (5.17)  and  the  definition  of  the 
advection  operator  (4.3),  Equations  (4.9)  and  (4.11)  can  directly 
be  translated  into  their  discretized  approximations: 


I  W  ^V,(r)  = 

Lc  on  3t  2  n  Lc 
r  sS  r  eS 

n  n 


-WUVV0(r  )-WV:UV„(r  )  + vWA  V. (r  ) 
on  2  n  on  2  n  on  2  n 


-  Y  WV  V_(r  )V0(r  )  -  Y  (WV  +WV  )Vi  (r  )V_(r  ) 
L c  onm  2  n  2  m  L  „  onm  omnv  I  n  2  m 
r  eS  r  eS 

m  m 


l  W  V  (r  )  =  l 

Ln  on  dt  3  n 
r  eS  r  eS 

n  n 


r  e  S 
o 


TJV-*  V  *  U-*  A 

-W  V_ (r  )-W  '  v _ (r  )  + vW  V7(r  ) 
on  3  n  on  3  n  on  3  n 


(5.18) 


•  Y  W  V7(r  ) V  (r  ) 

onm  3  n  3m 

r  eS 
m 

•  Y  (Wv  +Wv  ) (V  (r  ) +V  (r  ))V_(r  ) 
Lc  onm  omn  I  n  2  n  3  m 

r  eS 
m 

r  e  S 
o 


(5.19) 


Finally,  the  coefficient  matrix  on  the  left  hand  side  of  Equations 
(5.14),  (5.18),  (5.19)  is  inverted  over  to  the  right  hand  side, 
as  is  shown  here  for  Equation  (5.14): 


-rV.(r)  =  l  W"1 
3t  I  n  L„  on 

r  eS 
n 


tl.UV  Tl.v :  U  ,A 
-W  -  W  +  vW 
on  on  on 


r  eS 
m 


I  WV  V.  (r  ) 
L „  onm  1  m 


r0  e  S, 


(5.20) 
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with  W  denoting  the  inverse  from  the  left-hand  side, 
on 

This  system  (5.20)  (relating  to  Equation  (4.7))  is  the  final 
result  of  the  discretization  procedure.  Once  the  mean  flow  U  and 
appropriate  basis  functions  W  (y)  are  selected  and  thereby  the 
coefficient  matrices  determined,  it  can  be  analyzed  for 
bifurcations  and  numerical  solutions.  Therefore,  the  construction  of 
basis  functions  determining  incompressibility  and  boundary  conditions 
follows  next. 

6.  DIVERGENGE-FREE  BASIS  FUNCTIONS 

Using  Equation  (5.8),  a  velocity  field  V(r)  over  D  approximates 

like: 


V(r)  »  l  Wr  (r)V(r  ) 

r  eS  n 
n 


■  Vr) 


(6.1) 


The  basis  functions  W  (r) ;  r  e  S;  r  e  D  have  to  have  the  properties 

n 

of  weight  functions,  i.e.: 


W 


r 

n 


sup  W 
reD  Tn 


(r) 


1 


(6.2) 


W 

r 


n 


=  0; 


r  ,  r  e  S:  m  n 

n  m 


(6.3) 


The  condition  (6.3)  suggests  restricting  the  support  of  W  (r) 

n 

within  a  neighbourhood  of  r  . 

Although  a  rectangular  discretization  is  in  principle  possible 
and  seems  to  be  simpler  to  implement  (since  the  r^  e  S  would  sit  on 
a  rectangular  grid  naturally  induced  by  the  cartesian  coordinate  system 
as  well  as  rectangular  areas  of  supports  for  the  Wr^(r)),  it  leads  to 
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complications  with  the  overlapping  properties  of  the -supports  and 

condition  (6.3).  Most  applications  of  finite  element  techniques 

therefore  discretize  the  fluid  domain  into  simplices  (i.e.  triangles  in 

two  dimensions  and  tetrahedrons  in  three  dimensions)  with  the  set  S 

of  nodes  consisting  of  all  centroids  on  the  interfaces  between  the 

simplices.  The  area  of  support  for  W  (r)  is  in  this  case  the 

n 

union  of  the  two  simplices  sharing  the  interface  with  r  as  centroid. 

n 

See  Figure  11.  In  this  discretization,  the  individual  tetrahedron  is 
the  smallest  volume  element  permitted.  This  means  for  the 
incompressibility  condition  that  no  net  flux  is  permitted  through  the 
faces  F  enclosing  the  tetrahedron: 

/  V(r)n  dF  =  0  (6.4) 

F  h 


n  is  the  unit  vector  perpendicular  to  the  faces  F.  Using  Greens 

r 

Theorem  and  the  approximation  (6.1)  gives: 


/  V  (r)n„dF  =  /  div  V  (r) dr, dr0dr. 
'  n  F  '  n  l2: 


=  /  l 


T  i  8ri 


I  (r)V(rn) 


r  eS  n 
n 


drldr2dr3 


l  V  (r  )  / 

r  eS  T 

n 


3W  (r) 
r 

n 

3r . 


drldr2dr3 


(6.5) 


This  has  to  vanish  for  all  possible  V(r  ):  r  e  S.  Therefore: 

r  n  n 


/ 

T 


3W  (r) 
r 


n 


drldr2dr3 


=  0; 


dr. 

l 


i  =  1,2,3 


(6.6) 


jp 
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FIGURE  11 


Two  Tetrahedrons  Sharing  a  Face  with  Centroid  (Support 
for  Finite  Elements) 
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It  has  been  shown  by  Fortin  [11]  that  only  piecewise  continuous  basis 

functions  or  so-called  non-conforming  finite  elements  are  compatible 

with  the  incompressibility  condition.  The  simplest  possible  structure 

for  the  Wr  (r)  are  therefore  polynomials  defined  over  the  two  support 
n 

tetradehrons  individually  with  Wr  (r)  vanishing  over  the  rest  of  D. 

n 

That  means,  the  Wr  (t)  will  be  discontinuous  over  the  faces  of  their 

1  n 

support  tetrahedrons.  Conditions  (6.2),  (6.3),  (6.6)  determine  all- 

together  ten  conditions  per  tetrahedron  and  polynomial  since  the 

condition  W_  (r  )  =  sup  W„  (r)  can  be  written  as: 
rn  n  tl  r 

reD  n 


grad  W  (r  )  =  0 
r  n 
n 


or: 


3W  (r  ) 
r  v  n 
n 

3r. 

l 


=  0,  i  =  1 , 2 ,3 


(6.7) 


This  suggests  the  following  polynomial  for  W  (r) 

rn 


wr  (t) 
n 


C  +  Y  C.r.  +  V  C. .r.r . 
F  ii  . L .  liii 

~\  T  ^  J  J 


1<J 


(6.8) 


Now  a  linear  algebraic  system  of  ten  equations  can  be  set  up  to  solve 
for  the  ten  coefficients  C,  C^,  of  Equation  (6.8).  This  can  be 

considerably  simplified  by  the  introduction  of  barycentric  coordinates 
(see  Appendix  B)  \  ;  L  =  1,2, 3, 4  : 

Li 


r  •  t  A. 

iL  L 


r . 
i 


(6.9) 
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4 

1  =  l  At  (6.10) 

L  L 

Here  is  the  i-component  of  the  L-th  vertex  in  cartesian 

coordinates.  A  derivation  of  this  10  x  10  system  can  be  found  in 
the  Appendix  B. 

7.  DIRICHLET  BOUNDARY  CONDITIONS  AND  STRUCTURE  OF  MEAN  FLOW  U 

As  described  in  the  Appendix  B,  the  discretization  program 
partitions  a  given  rectangular  domain  D  in  fluid  space  into  rectangular 
elements,  each  one  of  which  is  subdivided  into  five  tetrahedrons  (see 
Figures  12a, b) .  Omitting  all  nodal  points  on  the  surface  planes  3D  of 
D  from  the  discretization,  i.e.,  setting  their  basis  functions  or 
finite  elements  identical  to  zero: 

W  (r)  =  0  (7.1) 

3D 

means  Dirichlet  boundary  conditions  for  the  perturbation  approximations 

V  : 
n 

Vn(r3D,t)  =  0  Yt,  (7.2) 

as  follows  directly  from  the  approximation  (6.1).  This  is  equivalent 
to  setting  the  V(r  )  in  (6.1)  equal  to  zero.  It  is  evident  that, 

o  U 

due  to  the  construction  of  the  discretized  ODE-systems,  these  Dirichlet 

-L 

boundary  conditions  hold  for  the  perturbations  ,  V.,  only  (i.e., 

to  their  discretization  respectively) ,  and  not  for  the  main  flow 

_N 

U,  which  can  go  through  D  or  3D  respectively  unaffected.  This 
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FIGURE  12A.  A-Configuration  of  Discretization 
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MGURE  12B.  B-Configuration  of  Discretization 
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creates  the  following  physical  situation  (see  Figure  13) :  Within  a 
stationary  flow  field  U  ,  say,  which  can  be  a  solution  of 

GA  L 

Equation  (4.1)  or  (4.4),  respectively,  there  exists  a  rectangular  domain 
D  within  which  U  can  be  approximated  by  U  .  The  stability 

GAL 

properties  of  U  within  D,  i.e.  the  development  of  a  perturbation 

_ V  _ ^ 

on  top  of  U,  is  determined  by  Equation  (4.7)  or  its  discretization 
(5.20)  respectively.  However,  due  to  the  Dirichlet  boundary  conditions, 
this  perturbation  V  cannot  propagate  out  of  D. 

For  U  ,  the  following  structure  is  built  into  the  discretization 
program: * 


E11  E12  E13 

U1  *  C1  +  CllV  *  C12X2  +  C13X3 


p  p  p 

21  22  23 

U2  =  C2  +  C21X1  +  C22X2  +  C23X3 


^31  ^32  ^33 

U3  =  C3  +  C31X1  +  C32X2  +  C33X3 


(7.3) 


For  the  exponents,  integers  between  1  and  3  inclusive  have  to  be 
chosen  as  input  for  the  discretization  program.  The  coefficients, 
however,  are  real  and  variable  parameters  of  the  system  (14.4).  Besides 
the  viscosity  v,  they  can  therefore  be  used  as  stability-determining 
parameters  for  bifurcations  or  as  functions  thereof. 


In  principle,  any  conceivable  structure  of  U  can  be  built  into  a 
discretization  program  for  (5.20).  The  structure  (7.3)  however  is 
relatively  easy  to  program  and  offers  a  considerable  choice  for 
bifurcation  parameters.  See  Appendix  B  for  details. 
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U 


ext 


FIGURE  13.  Flow  Through  Domain  D  by  U 


PART  3 


8.  INVESTIGATION  OF  BIFURCATIONAL  PROPERTIES 

To  test  the  system  (-5.14)  for  Hopf  bifurcations,  BIF0R2,  a 
FORTRAN  program  written  by  B.  Hassard*  is  used.  Since  the  CPU  time 

3 

of  the  program  is  roughly  proportional  to  N  ,  with  N  as  the  number 
of  dimensions  of  the  ODE-system,  the  lowest  possible  discretization 
had  to  be  used.  That  is,  a  domain  D  consisting  of  one  cube  with 
corner  coordinates  (0,0,0),  (0,0,1),  (0,1,0),  (1,0,0),  (0,1,1), 

(1,0,1),  (1,1,0),  (1,1,1)  is  split  into  five  tetrahedrons  according 
to  Figure  12  (A-configuration) .  Due  to  Dirichlet  boundary  conditions, 
only  the  4  nodes  on  the  interior  tetrahedral  faces  contribute  to  the 
system  (5.14),  which  therefore  becomes  12-dimensional. 

Since  BIFOR  2  is  explained  very  extensively  and  sufficiently 
in  theory  and  function  by  its  author  in  [5] ,  only  a  brief  introduction 
and  description  is  given  here. 

BIFOR  2  varies  the  bifurcation  parameter  n  in  the  Jacobian 
matrix  A(X*(n),ri)  of  the  ODE  system  until  the  set  of  eigenvalues 
Y(n),  i  =  1,...,N,  consists  of  one  complex-conjugate  pair  with  zero 
real  parts,  and  the  real  parts  of  the  rest  of  the  set  are  negative.  This 
is  performed  by  a  secant  iteration.  In  between  each  two  of  these 
secant  iteration  steps,  the  stationary  solution  X*(n)  for  the  new 
iterate  of  n  is  determined  by  Newton  iteration  before  the  new  set  of 
eigenvalues  is  evaluated.  Since  the  stationary  solution  of  the  system 
(5.14)  of  interest  is  always  the  trivial  solution  (n)  =  0,  Yn  the 

ie 
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user-supplied  initial  guess  for  the  stationary  solution  at  the  critical 
value  nc  of  the  bifurcation  parameter  is  already  known  to  be 
V  (y  )  =  0.  Once  X^(n  )  and  n.c  are  found,  the  program  calculates 
all  parameters  to  be  used  in  the  following  three  equations: 


X(t,n)  =  X*(y  )  + 


rrn, 


ReL2Tit/T(n)v 


ec 


T  (n)  =  — 

P 


1  +  T, 


— 1 

ZJ 

1 

zs 

1 

c 

.  y2  . 

(8.1) 


(8.2) 


8(n)  =  6. 


n-n. 


u. 


(8.3) 


These  are  first-order  approximations  for  the  periodic  perturbation 
X(t,n),  its  period  T(n),  and  the  Floquet  exponent  3(y),  which 
determines  the  stability  of  X(t,y): 


3  (n )  >  0  :  X(t,n)  is  unstable 
3(n)  <  0  :  X(t,n)  is  stable 


(8.4) 


y9  is  a  coefficient  from  a  Taylor  expansion.  It  determines  the 
criticality  of  the  bifurcation: 


>  0  :  supercritical  bifurcation;  X(t,n)  exists  on  RHS  of  y^ 
y^  <  0  :  subcritical  bifurcation;  X(t,n)  exists  on  LHS  of  y^ 

V  is  the  complex  eigenvector  of  the  Jacobian  matrix  A(X*(y  ) ,y  ) 
and  is  the  imaginary  part  of  the  complex  conjugate  pair  of 

eigenvalues  with  zero  real  part  at  y  =  y^  .  and  are 

coefficients  from  Taylor  expansions;  they  and  y9 


are  determined  as 
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analytic  functions  of  the  first  coefficient  of  the  Poincard  normal 
form  of  the  ODE  system. 


L/2 

5  =  X  Cn)  ?  +  l  C  ,(rO;|s 
j  =  l  3 


2j 


-  0  C I  5 


(€>n) 


(8.5) 


This  form  is  obtained  by  a  coordinate  transformation  of  the  system 
and  x  ,  B2j  U2  are  §iven  by: 


-Re(C  (0)) 

y2  =  a' (n  )  ’  e2  =  2Re Cc! CO) ) 

v  c 


(8.6) 


-(Im)  (C^O))  +u2m'(Tic)) 


(8.7) 


with  the  first  coefficient  of  (8) : 


CJ°>  =  if 


o 


(2Ulf20) f2Ulfll)  '  2 


T  i  2  1  i  T-  ,2 

2U  f  -  —  2U  f 
111 1  3  1  1  201 


21 


(8.8) 


here  the  are  transposes  of  the  left  eigenvector  of  the  Jacobian 

matrix  A(X*  (n£)  ,  nc)  for  the  eigenvalue  X^Oi  )  =  and  f?Q,  f^, 

G  are  partial  derivatives  of  the  right-hand  side  of  the  ODE  system. 
Moreover: 


X{(nc)  =  a'Cnc)  +  im'(nc)  (8.9) 

is  the  eigenvalue  X^(n)  derived  after  q  at  q  .  The  program 
evaluates  these  equations*  and  thereby  supplies  the  numerical  values 


*Except  for  the  Poincare  normal  form  (8.5),  since  only  its  coefficient 
C  (0)  is  used. 


65 


for  all  the  parameters  in  Equations  (8.1),  (8.2),  (8.3).  Their  complete 
derivation  can  be  found  in  [5] . 

The  first  flow  to  be  tested  was: 


U  =  0;  U  =  C  X;  U=0 
x  y  yx  z 


(8.10) 


with  :  n  =  C 

yx 

Around  80  runs  were  performed  and  for  bifurcation  parameters  the 

coefficient  C  or  C  were  used.  After  initial  trial  runs,  it 
x  yx 

turned  out  that  the  flow  was  unstable  for  positive  viscosities 

throughout.  That  is,  BIFOR 2  could  not  zero  any  of  the  real  parts 

of  the  eigenvalues;  they  were  all  positive  and  bounded  away  from  zero. 

Since  this  behaviour  was  observed  for  all  flows  tested  so  far, 

negative  viscosities  had  to  be  used  throughout.  For  the  flow  (8.10) 

viscosity  values  between  -.01  and  -5  were  used.  For  a  viscosity 

of  -1,  the  flow  became  unstable  at  C  ^6.56.  However,  the 

yx 

_ v 

bifurcation  was  not  Hopf  and  the  perturbation  therefore  stationary. 

This  was  found  for  all  bifurcations  for  the  flow  (8.10).  Moreover, 
it  turned  out  for  this  and  all  other  tested  flows  that  the  critical 
value  for  the  bifurcation  parameter  is  proportional  to  the  absolute 
value  of  the  viscosity: 

nc  ^  | v |  (8.11) 

The  phenomenon  of  stationary  bifurcation  shows  that  the  flow  (8.10) 
somehow  changes  its  structure  before  it  can  undergo  a  Hopf  bifurcation. 
Therefore,  nonzero  initial  guesses  for  the  perturbation  V 


were 
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experimented  with  to  find  a  Hopf  bifurcation  for  a  flow  of  the  type 

(8.10)  +  (n  )  .  However,  the  program  always  either  converged  to 

— 

the  trivial  solution  e  0  or  error-terminated  in  the  Newton 
iterations.  Moreover,  the  additional  Newton  iterations  increased  the 
CPU  time  drastically  and  this  approach  was  therefore  abandoned  for 
all  other  tested  flows  as  well. 

A  listing  of  all  the  other  tested  flows  with  a  discussion  of 
their  bifurcational  properties  is  given  in  Appendix  C.  Here,  only 
the  main  characteristics  are  summarized: 

1.  As  mentioned  already,  only  negative  viscosities  can  be  used.  All 
tested  flows  seem  to  be  unstable  for  positive  viscosities.* 

2.  Flows  with  linear  or  no  shear  or  one  component  only,  have  only  real 
eigenvalues.  As  a  rule,  the  program  can  find  the  critical  value 

of  the  bifurcation  parameter,  where  the  largest  eigenvalue  vanishes. 
This  can  be  understood  as  an  inability  of  no  shear  or  linear  shear 
to  develop  periodic  perturbations.  These  types  of  profiles  have  to 
reset  themselves  (e.g.,  by  a  stationary  bifurcation)  to  nonlinear 
ones  before  they  can  undergo  a  Hopf  bifurcation. 

3.  The  critical  value  nc  and  the  frequency  u)Q  are  proportional  to 
the  magnitude  of  the  viscosity,  whereas  3^  and  y9  are  inversely 
proportional . 

4.  Most  of  the  flows  show  qualitatively-  correct  behaviour;  an  increase 
of  shear  in  the  mean  flow  U  leads  to  its  instability.  In  the 
case  of  linear  shear  this  occurs  by  a  stationary  bifurcation  and 


*See  Appendix  C  for  an  exception. 
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in  the  case  of  nonlinear  shear  a  supercritical  stable  Hopf 
bifurcation  occurs. 

To  obtain  a  superposition  of  three  frequencies  leading  to  a 
3-torus  as  the  invariant  manifold  in  solution  space,  the  following 
idea  was  applied  first. 

The  analytic  approximation  for  V^(r,t)  according  to 
Equation  (8.1)  is  evaluated  for  a  set  of  fixed  time  points  t^, 
i  =  1,...,N,  say,  and  the  obtained  values  V  (r,t^)  are  substituted 
into  the  system  for  V?(r,t)  (i.e..  Equation  (5.18))  and  tested  for 

bifurcations  one  after  the  other.  This  way  one  obtains  N  secondary 
orbits  for  the  second  periodic  perturbation  V  ^(r,t),  i  =  1,...,N, 
again  by  Equation  (8.1).  That  is,  the  2-torus  relating  to  two 
bifurcations  is  discretized  into  N  closed  orbits.  Substituting  the 
V^(r,t^)  and  evaluations  at  fixed  time  points  for  the  secondary 
orbits  V  . (r,t.),  j  =  1,...,M,  say,  into  the  system  for  V  (r,t) 

(i.e..  Equation  (5.19))  and  testing  for  bifurcations  one  obtains  N  x  M 
orbits  for  the  third  bifurcation,  i.e.,  each  secondary  orbit  is 
discretized  into  M  points  which  degenerate  into  orbits  for  the  third 
bifurcation.  This  N  x  M  discretization  of  the  3-torus  can  now  be 
used  to  approximate  any  trajectory  on  it  by  use  of  interpolation 
methods.  However,  this  approach  was  abandoned  since  for  the  flows 
tested  secondary  and  third  order  bifurcations  could  not  be  obtained 
for  each  different  time  setting  on  the  first  orbit  and  secondary  orbits. 
Moreover,  it  is  not  established  yet  whether  this  type  of  discretization 
in  time  is  legitimate  for  the  testing  of  non-autonomous  systems  for 
higher-order  bifurcations.  On  the  other  hand,  the  flows  which  showed 
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such  a  sequence  of  a  first,  secondary  and  third  order  bifurcation  at 
least  for  some  time  settings  turned  out  to  have  better  recurrency 
properties,  as  was  found  out  by  the  second  approach,  which  combines 
BIFOR 2  with  ODE  solvers. 

9.  APPLICATION  OF  ODE  SOLVERS 

Since  the  construction  of  a  discretized  torus  cannot  be  pursued 
because  of  the  reasons  mentioned,  numerical  solutions  generated  by  ODE 
solvers  offer  a  reasonable  alternative.  For  this  method,  the  non- 
autonomous  system  for  V^Crjt)  was  used  with  the  approximation  by 
Equation  (8.1)  for  V^(r,t)  determined  by  BIFOR 2  . 

Both  predictor-corrector*  and  Runge-Kutta**  solvers  were  applied. 
Unfortunately,  the  corrector-predictor  method  had  to  be  abandoned 
despite  the  favourable  CPU  times  (which  are  about  one  fourth  of 
the  Runge-Kutta  programs)  since  it  caused  in  all  its  options  and 
variations  a  '’sawtoothing"  on  the  solutions,  as  can  be  seen  by 

Figures  Tl,  T2  in  Appendix  F.  However,  the  CPU  time  on  the  single¬ 
precision  RKF45  routine  could  be  dropped  to  about  one  sixth  of  the 
double-precision  DVERK  routing  by  recompiling  the  user-supplied 


IMSLDPLIB  DGEAR:  Gears  or  Adams  corrector-predictor. 

IMSLDPLIB  DVERK:  Sixth  order  Runge-Kutta-Verner  method. 

RKF45:  Fourth-fifth  order  Runge-Kutta-Fehlberg  method. 

Authors:  L.F.  Shampine,  H.A.  Watts;  Sandia  Corporation, 

Albuguergue,  New  Mexico  87115. 

See  [29]  for  a  description  of  RKF45. 
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routine  which  evaluates  the  ODE- system  on  a  FORTRAN  77  compiler  with 
highest  efficiency  option.* 

The  following  four  flows  were  used  for  the  described 
application  of  ODE  solvers.: 


-1. 


(I) 


This  flow  bifurcates  stable  and  subcritically  at  n  =  .872  with  a 


period  of  T^  =  3.24.  Higher  order  bifurcations  using  the  described 
method  of  torus  discretization  could  not  be  found  except  for  a  station¬ 
ary  one  with  V1(r,t  =  0)  and  at  n  «.952.  Numerical  solutions  were 
attempted  at  n  =  .9,  .915,  .93,  1  with  initial  conditions  equal  to 
.001  for  all  twelve  velocity  components  of  V2(r,t=0).  For  n  =  .9 
the  flow  was  completely  stable  in  its  first  bifurcation.  The  running 
time  was  T  =  30  (300  time  steps  with  .1  for  step  size)  and 

V- 

vanished  completely.  For  n  =  .915,  however,  V?(r,t)  grew 
unbounded  with  little  sign  of  recurrency  and  therefore  error- 
terminated  at  T  =  6.2.  Essentially  the  same  behaviour  was  observed 
for  n  =  .93  and  1  with  an  error- termination  at  T  =  5.3  and  2.6  . 


U  =  5nY2;  U  =  3nZ2;  U  =  nX2 
x  y  z 


(II) 


v  =  -1  . 


CPU  times  on  the  U/A  AMDAHL  were  by  order  of  magnitude  ca.  150 
seconds  for  1000  timesteps  with  stepsizes  of  ca.  .02-  .05 
simulation  time  on  DVERK  IMSLDPLIB  for  both  random  perturbed  and 
unperturbed  runs.  This  could  be  reduced  to  ca.  60  sec.  for  random 
perturbed  runs  and  ca.  25-  30  sec.  for  unperturbed  runs  on  RKF45. 


*  * 


All  times  T  or  T  are  understood  as  simulation  time. 

P 
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The  first  bifurcation  occurs  (stable  and  supercritical)  at  n  =  1.488 
with  a  period  of  =  .23  .  A  secondary  Hopf  bifurcation  occurs  at 
n  =  1.881  for  V^(r,t  =  0).  ODE  runs  were  tried  at  n  =  1.9,  2.23, 
2.25,  2.3,  2.4  with  V  (r,t = 0)  =  .0001  for  all  twelve  components. 

For  n  =  1.9,  the  primary  orbit  was  still  stable  despite  the  secondary 
bifurcation  at  n  =  1.881;  after  T  =  30,  the  twelve  components  of  V 
were  of  an  order  of  magnitude  of  around  10  ^ .  For  n  =  2.23,  2.25, 
2.3,  2.4,  the  run  error-terminated  at  T  =  14.42,  21.25,  11.19,  4.6 
due  to  unbounded  growth  of  V  .  For  none  of  these  parameter  settings 
did  V9  show  any  significant  periodicity  or  recurrency,  and  the 
unbounded  growth  causing  the  error  termination  occured  rather  suddenly. 

(Ill)  U  =  -9n  +  (6  +  9n)Y2 

X 

=  -6n  +  (5  +  6n)Z2 

u  =  -n  +  (3  +  n)x2 
z 

v  =  -1. 

This  flow  undergoes  the  first  Hopf  bifurcation  at  n  =  .189,  the 

V 

secondary  one  for  V2(r,t  =  0)  at  n  =  .637  and  the  third  order  one 
for  V^(r,t  =  0)  at  n  =  1.003.  The  affine  functional  dependencies 
of  the  shear  coefficients  on  the  bifurcation  parameter  r\  seem  to 
increase  the  chances  for  bifurcations  up  to  third  order  considerably. 
Amongst  the  tested  flows  with  pure  linear  dependencies,  none  could  be 
found  with  bifurcations  up  to  third  order. 

The  tested  parameter  values  for  this  flow  were  n  =  .4,  .5,  .52, 
.53,  .55,  .6,  .7  .  The  Figures  Tl,  T2  (Appendix  F)  show  projected 
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trajectories  for  +  V0  at  n  =  .55  and  a  period  of  T^  =  1.66  for 
V  .  Strong  recurrency  is  evident,  but,  as  can  be  concluded  by  the 
graphs,  unbounded  growth  occurred  which  error-terminated  the  run  at 
T  =  24.523.  The  trajectories  for  all  other  parameter  settings  were 
similar  with  unbounded  growth  and  error- termination  occurring  for 
n  =  .6  and  .7  at  T  =  9.4  and  7.44  respectively.  For  the  settings 
of  n  lower  than  .55  unbounded  growth  did  not  occur  within  the 
chosen  time  settings,  which  went  as  high  as  T  =  40  for  n  =  .52. 

(IV)  U  =  -9 a  +  (4  +  9n)Y2 

X 

=  -6n.  +  (5  +  6n)Z2 
U  =  -n  +  (3  +  n)X2 

Zj 

V  =  -1  . 

Although  the  only  difference  between  this  flow  and  flow  (III)  is  in 
the  affinity  constant  of  the  shear  coefficient  in  the  U  -component, 

X 

the  trajectories  of  the  solution  curve  looked  somewhat  more  promising  as 
far  as  the  resemblance  to  projections  of  2-  and  3-  tori  is  concerned. 
Therefore,  it  has  been  tested  much  more  extensively.  The  flow 
bifurcates  at  n  =  .242  and  at  n  =  .329  and  .347  for  V  (r,t=0) 
and  V2(r,t  =  0)  respectively.  ODE  runs  were  performed  between 
n  =  .28  and  .5  .  All  graphs  (Appendix  F)  show  projections  of  the 
+  V0  solution  curve  onto  the  three  velocity  component  planes  for 
each  one  of  the  four  monitored  points  in  D  . 

Following  these  graphs  from  lower  to  higher  values  of  n  with 
initial  conditions  for  V0  close  to  zero,  it  can  be  seen  that  the 
orbit  for  is  being  followed  closely  (i.e.,  stays  small)  for 
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a  number  of  revolutions  which  decreases  with  increasing  n  ;  e.g.,  for 
D  =  .3  (Figure  T3)  this  primary  orbit  is  being  closely  retraced  for 
ca.  6  times,  whereas  at  n  .45  (Figure  T4)  only  for  ca.  1^  times. 
After  these  initial  revolutions  a  number  of  transitional  revolutions 
follows  (again  decreasing  with  increasing  p),  where  the  system  passes 
from  the  V^-orbit  to  its  new  invariant  manifold.  For  p  =  .3 
(Figure  T3) ,  this  invariant  manifold  might  not  have  been  reached  within 
a  time  T  =  30.  For  all  settings  of  .325  <_  p  <  .449,  (Figures  T5, 

T6,  T7,  T8)  the  running  time  was  sufficient  to  reveal  a  limit  cycle 
as  invariant  manifold;  except  for  p  >_  .449  (Figure  T4)  where  the 
solutions  start  to  grow  unbounded.  The  shape  of  a  double  loop  of 
this  limit  cycle,  which  can  be  seen  by  most  of  these  projections, 
suggests  period  doubling.  It  could  however  by  determined  by  setting 
marks  on  the  trajectories  after  every  period  T  (p)  as  calculated 
by  Equation  (8.2)  that  this  period  stays  roughly  the  same  for  + V0  . 

For  increasing  p  the  limit  cycle  becomes  more  attracting.  It 
changes  and  grows  in  size  up  to  p-  .449  (Figure  T4) ,  where  it  seems  to 
"fan  out"  in  a  particular  section  due  to  an  increase  of  the  X-component 
on  top  of  it's  periodicity,  as  can  be  seen  by  projections  into  the 
XY-  and  XZ- planes.  The  running  time  of  T  =  .32  however  was  not 
sufficient  to  determine  whether  the  X-components  will  grow  unbounded 
or  approach  a  finite  limit.  Increasing  p  by  .001  to  p  =  .45* 
causes  an  error-termination  due  to  unbounded  growth  of  the  X-components 
at  T  =  22.82.  Moreover  a  spreading  of  the  trajectories  can  be 
observed  not  only  in  the  sections  of  this  error  termination,  but  in 
other  sections  as  well.  To  observe  some  of  the  behaviour  of  the  flow 

•k 

See  Figure  T4  (p  =  .449),  which  looks  identical  (except  for  unbounded 
growth) . 
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for  p  >_  .45  (Figure  T9),  after  an  error  termination  due  to  unbounded 

growth,  the  last  determined  value  of  V  was  reduced  by  multiplying 

-4  -6 

it  with  a  factor  in  the  orders  of  magnitude  of  10  to  10  to 
continue  the  run.  Evidently  the  legitimacy  of  this  is  questionable, 
since  beyond  p  =  .45,  the  flow  seems  to  be  in  a  region  where  it 
cannot  be  described  any  more  by  the  present  model.  The  assumption  for 
this  reduction  trick  is  that  the  model  still  works  reasonably  well  close 
to  the  time  point  of  sudden  unbounded  growth  even  for  p  >_  .45  in 
case  of  the  present  flow,  and  it  can  therefore  be  understood  as  a 
sequence  of  different  runs  with  different  initial  conditions.  Therefore, 
the  graphs  for  p  =  .47  (Figure  T9)  show  wider  spacing  of  trajectories 
and  a  certain  randomizing  effect  due  to  this  reduction  method  is 
evident  too.  Moreoever,  Y-  and  Z-components  can  grow  unbounded  as 
well  for  these  higher  parameter  settings,  as  the  YZ-pro j ection  in 
particular  shows. 

In  case  of  the  present  "one  cube  discretization"  the  coordinates 
of  the  four  points  in  fluid  space  are: 

Point  1  (1/3,  1/3,  1/3) 

Point  2  (1/3,  2/3,  2/3) 

Point  3  (2/3,  1/3,  2/3) 

Point  4  (2/3,  2/3,  1/3) 

according  to  the  A-configuration  (see  Appendix) .  By  the  amplitude 
(i.e.,  "size")  and  shape  of  the  trajectories  of  their  associated 
velocity  vectors,  these  four  points  fall  into  two  groups.  Group  one, 
consisting  of  point  1  and  3  has  orbits  which  are  ca.  one  to 

two  orders  of  magnitude  larger  than  those  of  group  two,  i.e.,  point 
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2  and  4.  This  seems  to  be  the  reason  why  the  limit  cycles  of 
V^(r= Point  l,t)  +  V =  Point  l,t)  and  V^ (r = Point  3,t)  + 

V?(r=  Point  3,t)  stay  within  the  order  of  magnitude  of  the  orbits  of 

—2k  _2». 

(r = Point  l,t)  and  (r = Point  3,t)  respectively  up  to  the 
occurrence  of  global  instability  at  n  =  .449  (Figure  T4)  whereas  the 
ratio  of  these  respective  orbits  and  limit  cycles  for  group  two  is  at 
least  one  order  of  magnitude.  For  group  two  all  six  projections  show 
the  double  loop  in  the  limit  cycle  whereas  for  group  one  it  can  be 
seen  only  in  the  XY  projection  of  point  one  and  the  YZ  projection 
of  point  three. 

10.  SIMULATION  OF  RANDOM  PERTURBATIONS  AND  THEIR  EFFECT  (ON  FLOW  IV) 

To  simulate  the  aforementioned  "appropriate"  combination 
of  perturbations  in  connection  with  the  use  of  ODE-solvers,  IMSLLIB 
random  number  generators  are  being  used  on  Flow  (IV) .  The  four  points 
are  being  perturbed  independently  of  each  other  in  the  following 
way:  Routine  GGSPH  supplies  the  coordinates  of  a  point  on  the  unit 

sphere  random  determined  under  the  assumption  of  uniform  probability 
distribution  over  the  sphere.  These  coordinates  are  being  multiplied 
with  one  and  the  same  factor,  which  is  random  generated  by  the  routine 
GGNQF,  a  generator  of  Gaussian  distributed  variates  with  mean  zero  and 
a  user-determined  variance  V  .  The  three  component  perturbations 
obtained  this  way  are  being  added  to  the  appropriate  velocity 
components  of  the  point  considered  perturbed.  After  that,  GGNQF  is 
used  again  to  random-generate  the  number  of  unperturbed  timesteps  to 
follow.  The  average  of  this  number,  i.e.,  the  variance  V  ,  is  again, 
as  for  V  ,  set  by  the  user.  This  time,  only  the  absolute  integer  part 
of  the  generated  number  is  used,  since  timesteps  are  positive  and 


75 


integer  only.  That  is,  the  number  of  perturbations  decreases  with  in¬ 
crease  of  the  temporal  variance.  The  main  characteristic  of  the  system 
under  random  perturbations  is  a  global  instability  for  all  tested  bi¬ 
furcation  parameter  values.  For  a  setting  of  n  =  .32  (Figure  Til) 
and  a  spatial  variance  of  as  low  as  =  .005  with  a  temporal  variance 

of  V  =  2  ,  most  of  the  projections  still  showed  a  strong 
resemblance  to  the  ones  without  random  perturbations  at  n  =  .32 
(Figure  T10) .  The  run,  however,  error-terminated  due  to  unbounded 
growth  at  T  =  11.35.  On  the  other  hand,  a  run  of  the  same  bifurcation 
parameter  setting  and  temporal  variance  and  =  .03  showed  a  rather 
high  degree  of  randomness  especially  for  the  projections  of  the 
points  of  group  two  but  would  not  error-terminate  within  the  total 
running  time  of  T  =  15.  However,  on  the  basis  of  the  runs  performed 
so  far  probability  assumptions  about  stability  of  a  certain  combination 
of  bifurcation  parameter  and  random  variances  still  seem  to  be 
reasonable.  That  is  for  lower  n  and  and  higher  V  ,  the  chances 

for  a  "more  stable"  behaviour  are  better.  To  overcome  this  problem 
of  global  instability,  the  reduction  method  as  described  earlier  was 
again  employed.  (In  general,  the  projections  of  group  two  show  a 
higher  degree  of  randomness  than  group  one.) 

Although  the  spatial  structure  of  the  random-perturbed 
trajectories  is  in  most  cases  quite  similar  to  the  unperturbed 
trajectories,  their  temporal  structure  is  quite  different  and  shows 
more  randomness.  For  example,  there  are  regions  where  sets  of 
trajectory  segments  are  organized  so  that  they  appear  nearly  parallel 
to  each  other.  However,  these  sections  are  not  being  traced  in  some 
consecutive  order,  but  rather  randomly,  as  can  be  seen  when  the 
trajectory  is  being  drawn  on  a  graphics  device. 
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11.  SPECTRAL  ANALYSIS  (OF  FLOW  IV) 

To  obtain  a  comparison  with  the  inertial  subrange  law  (Appendix  E), 
fast  Fourier  transforms  were  performed  both  on  most  of  the  random- 
perturbed  ODE-runs  as  well  as  on  the  unperturbed  ODE-runs. 

First,  the  magnitudes  V(r^,t^)  of  the  total  velocity  vectors 
for  each  timestep  t;  and  each  of  the  four  monitored  points  r^  , 
j  =  1,...,4  are  determined  by  the  (obvious)  formula: 


V(r . ,  t . ) 

j  i 


k=l 


(vstat 


(11.1) 


with  k  counting  through  the  X, Y, Z-components  and  V  =  velocity 
from  steady  input  profile;  =  first  perturbation  from  bifurcation; 

V?  =  second  perturbation  growing  on  and  determined  by  ODE-solvers. 

After  that,  the  mean 


V(r.) 
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(11.2) 


taken  over  all  i  timesteps  is  subtracted  from  each  individual 

max 

|V(r^,t^)|,  to  obtain  a  pure  time- dependent  part  of  the  velocity: 


V  (r . ,  t .  )  I  e  V  (r  . ,  t .  )  =  I  V(r  .  ,  t . )  I  -  |V(r.) 

p  j  l  1  p  j  l  1  j  i  1  1  j 


(11.3) 


This  way,  a  possible  nonzero  mean  of  stemming  from  its  nonsymmetric 

trajectory  w.r.t.  zero  is  also  subtracted. 

The  velocities  are  now  analyzed  for  their  power  spectra  in 
two  pairs  of  two  points  (1,2  and  3,4)  by  the  IMSLLIB  routine  FTFPS 
to  obtain  frequency  power  spectra  and  cross  spectra  estimates.  For 
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this  spectral  analysis,  the  entire  time  series,  i.e.,  a  particular 

ODE-run  in  this  case,  is  subdivided  into  intervals  of  equal  size, 

which  must  be  a  power  of  two,  and  whose  number  must  divide  the  number 

of  samples  in  the  entire  time  series  evenly.  If  this  is  not  possible, 

the  complete  number  in  the  time  series  is  increased  by  samples  of 

value  zero,  commonly  referred  to  as  "padding".  For  example,  in  the 

case  of  1000  timesteps,  which  was  the  largest  and  most  frequently 

6  7 

used  number  of  timesteps,  subinterval  sizes  of  64  (=  2  )  or  128  (=  2  ) 
were  chosen.  A  time  series  of  size  1024,  contains  exactly  16  or  8 
respectively  of  these  subsamples.  A  run  of  1000  timesteps  was 
therefore  padded  with  an  extra  24  timesteps  of  zero  value. 

Each  one  of  these  subsamples  is  Fourier-analyzed  and  an  average 
spectrum  is  determined  from  all  these  subsample  spectra.  Clearly,  this 
causes  a  smoothing  in  this  average  spectrum,  which  is  necessary  to 
determine  a  definite  slope  or  trend,  that  would  be  obscured  by  dense 
oscillations  in  a  spectrum  obtained  without  averaging. 

To  obtain  a  direct  comparison  with  the  inertial  subrange  law 
(Appendix  E)  the  spectrum  frequencies  are  rescaled  as  wavenumbers  under 
the  assumption  of  the  Taylor  hypothesis.  The  diagrams  show  the  power 
spectra  (i.e.,  spectral  energy  density)  on  the  ordinate  versus  wave- 
numbers  on  the  abscissa,  with  both  axes  on  decadic  log  scale. 

A  spectral  analysis  was  performed  for  the  following  eleven 
parameter  settings:*  n  =  . 27(p)  (T12,S12),  .28(p)  (T13,S13),  .3(6) 
(T3,S3,T14,S14) ,  .325(b)  (T5,S5,T15,S15) ,  .33(p)  (T16,S16), 

•k 

See  Appendix  F  for  Figures;  T  (trajectories),  S  (spectra), 
p  =  random  perturbations  superimposed. 

U  =  without  random  perturbations, 
b  =  both  with  and  without  random  perturbations. 


* 


78 


. 35 (U)  (T6,S6),  • 4 (U)  (T7,S7),  .445(p)  (T17,S17),  .447(U)  (T8,S8), 

.449(11)  (T4,S4),  .47(U)  (T9,S9).  For  lower  parameter  settings,  such 
as  n  =  .27, .28,  the  spectra  show  an  almost  horizontal  section  for 
wavenumbers  around  .2  to  .5  and  spectral  densities  around  .3  for 
group  I  (points  1,3)  and  wavenumbers  around  .2  to  .7  with  spectral 
densities  around  .07  to  .001  for  group  2  (points  2,4).  This 
difference  of  ca.  one  order  of  magnitude  between  the  spectra  of  the 
two  groups  reflects  the  already  mentioned  difference  in  size  of  one  order 
of  magnitude  between  the  limit  cycles  of  the  two  groups.  This  section 
in  the  spectrum  can  clearly  be  associated  with  the  production  range  of 
large  eddies  in  a  turbulent  flow. 

After  this  initial  section,  the  spectra  fall  off  with  a  slope  of 

roughly  -7  to  wavenumbers  of  ca.  2  to  3  and  spectral  densities 
-4  -5 

around  10  to  10  .  The  ranges  between  this  section  and  wavenumbers 
of  ca.  8  to  9  shows  more  of  an  oscillatory  behaviour  in  their  spectra, 
which  is  due  to  the  logarithmic  scale  and  therefore  a  higher  density  of 
data  points  on  the  graph.  However,  the  mean  slope  in  these  sections 
compares  quite  well  to  -5/3. 

Increasing  n  to  .3  changes  the  slope  of  the  midsection  from 
ca.  -7  to  ca.  -5  and  a  slight  steepening  of  the  mean  slope  in  the 
right  section  for  large  wavenumbers  from  ca.  -5/3  to  ca.  -5/2  . 

This  spectrum  stays  essentially  the  same  for  n  =  .325  .  Also,  the 
unperturbed  spectra  for  n  =  .3,  .325  show  in  principle  the  same  shape 
as  their  perturbed  counterparts  just  described. 

The  slopes  closest  to  -5/3  can  be  found  for  the  spectra  of 
group  1  at  h  =  .33  and  .445.  An  analysis  of  these  spectral  slopes 
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is  displayed  in  Table  1  together  with  the  ones  of  Group  2. 

In  general.  Group  2  again  has  lower  spectral  densities  than  Group  1. 

The  concavity  of  Group  2  at  the  lower  spectral  end  indicates  an  eddy 
production  range.  Although  point  1  shows  an  overall  average  slope 
of  -1.7  for  n  =  .33  and  the  same  in  the  lower  wavenumber  range  for 
n  =  .445,  it  must  be  noted  that,  for  definite  answers  w.r.t.  acceptance 
of  these  results,  they  would  have  to  be  reproduced  with  sample  sizes  of 
at  least  an  order  of  magnitude  larger  than  the  presently  analyzed  ones. 
Nevertheless,  the  change  in  shape  and  magnitude  of  the  spectra  between 
n  =  .325  and  .33  is  remarkable,  especially  since  it  cannot  be 
concluded  from  the  trajectory  graphs.  The  Table  1  clearly  shows 
for  the  points  of  Group  1  spectral  slopes  in  general  closer  to  -5/3 
than  the  points  of  Group  2.  This  could  have  something  to  do  with  the 
Group  2  points  being  in  a  region  of  shear  high  enough  so  that 
generation  of  homogeneous  and  isotropic  turbulence  cannot  take  place. 

Looking  at  the  unperturbed  spectra  for  n  .35,  the  most  evident 
difference  w.r.t.  the  already  discussed  spectra  are  the  strong 
oscillations  spanning  ca.  2  orders  of  magnitude.  These  can  be 
expected  on  theoretical  grounds,  since  the  associated  trajectories  are 
periodic  and  therefore,  only  a  discrete  set  of  frequencies  or  wavenumbers 
respectively  can  occur.  The  spectra  for  n  =  .35  and  .4  look 
essentially  the  same.  The  drastic  decrease  of  oscillations  for  Group  1 
and  complete  disappearance  of  spectrum  for  Group  2  at  higher  wavenumbers 
(around  1  to  5)  indicates  the  absence  of  high  frequencies,  which  can  be 
concluded  from  the  associated  trajectory  graphs.  At  n  =  .447  and  .449 
the  spectra  compare  again  quite  well.  Compared  to  the  ones  at  n  =  .35 


TABLE  1 


ANALYSIS  OF  SPECTRAL  SLOPES  FOR  n  =  .33  (S16)  AND  .445  (S17) 


n 

.33 

.445 

Number  of  Timesteps; 

445,  .03 

500,  .03 

btepsize 

Mean  Magnitude  of 

.01 

.01 

Random  Perturbations  (V  ) 

s 

Mean  Number  of  Unperturbed 

2 

2 

Consecutive  Timesteps  (V  ) 

L 

Sample  Size/Subsample  Size 
=  Number  of  Subsamples 

448/64  =  7 

512/64  =  8 

-2.6C.3-  1.5); 

-1 . 7  ( . 3  -  1.5) ; 

Mean  Slopes 

Point  1 

(with  their 
respective 

-  .8  (1.5  -  10) 

-  1  (1.5-6) 

wavenumber 

ranges) 

Point  2 

.  5  (.3  -  1); 

-1.1(1  -  5) 

-1.3(1  -  5) 

-  2  (.2-1); 

-1 . 5 ( .  2  -  1); 

Point  3 

-  .7(1-4) 

-  .7(1  -6) 

-2 . 5 ( .  3  -  .6) 

Point  4 

(see  graph) 

-  2  (6-1) 
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and  .4,  the  oscillations  are  denser  and  span  a  wider  range  of 
spectral  density,  which  has  its  maximum  with  close  to  8  orders  of 
magnitude  for  points  3  at  n  =  .447  .  In  general,  a  shift  to 
higher  spectral  densities  can  be  observed,  which  is  a  natural  consequence 
of  the  increase  in  amplitudes  of  and  with  increasing  n •  The 

appearance  of  strong  spectral  oscillations  for  Point  1  in  the  wavenumber 
range  between  2  and  6  is  also  conspicuous,  since  it  cannot  be  inter¬ 
preted  from  any  particular  structure  in  the  trajectory  graphs. 

A  remarkable  flattening  of  mean  spectral  slope  and  decrease  in 
amplitudes  of  spectral  oscillations  can  be  seen  at  the  highest  parameter 
setting  analyzed,  i.e.  at  n  =  .47.  The  associated  trajectory  actually 
consists  of  10  short  trajectories  with  an  average  of  100  timesteps 
each.  Here,  the  already  discussed  reduction  technique  was  employed  to 
deal  with  the  unbounded  growth  of  .  This  could  explain  the  flattening 
of  spectral  slope,  which  covers  only  up  to  ca.  one  and  a  half  orders 
of  magnitude  in  spectral  density.  It  is  possible  that  in  this  range 
of  the  parameter  n  the  model  loses  its  applicability  due  to  un¬ 
realistically  high  shears  in  the  fluid,  which  cannot  be  modelled  on  the 
basis  of  NS-equations ,  or  at  least  not  by  the  finite  element  technique 
here  employed  and  rather  coarse  discretization. 

12.  SUMMARY  AND  CONCLUSIONS 

To  begin  with  the  theoretical  mechanism  of  the  model,  its  main 
characteristic  probably  lies  in  the  fact  that  the  NS-equations  can  be 
considered  as  being  both  valid  and  invalid  within  its  framework.  Locally, 
they  are  valid,  since  the  system  follows  sections  of  NS-solutions 
within  a  certain  neighbourhood  of  every  point  in  state  space. 
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Globally,  they  are  not  valid,  since  the  perturbed  trajectory  through 
state  space  of  the  system  cannot  be  considered  a  NS-solution,  if 
monitored  over  a  time  period  long  enough,  since  the  system  transfers 
between  sections,  which  are  not  connected.  The  trajectory  of  the  system 
is  therefore  not  identical  to  a  solution  curve  of  the  NS-equation.  This 
distinction  between  local  and  global  applicability  is  connected  to  a 
distinction  between  two  types  of  vectors  (or  spaces,  within  which  these 
vectors  are  defined)  upon  which  part  of  the  model  is  based. 

In  any  physical  context,  a  vector  has  three  cartesian  components 
in  its  general  form,  unless  the  coordinate  system  can  be  rotated  and 
these  three  components  are  reduced  in  number  in  adapting  to  the  partic¬ 
ular  geometry  at  hand.  The  most  important  characteristic  of  these  three 
components  is  that  they  are  of  the  same  type,  e.g.,  velocity,  accelera¬ 
tion  or  force.  Moreover,  they  are  all  associated  with  one  and  the  same 
point  in  physical  three-dimensional  space.  Vectors  of  this  type  could 
be  considered  ’’physical"  vectors,  and  they  are  local,  since  they  apply 
to  only  one  particular  point  of  a  system  under  consideration. 

On  the  other  hand,  any  point  in  any  axis  diagram  can  be  considered 
a  vector.  If  the  two  axes  on  the  diagram  represent  different  physical 
quantities,  any  rotation  (except  by  multiples  of  tt/2)  will  lead  to  a 
new  diagram  with  its  new  axes  as  linear  combinations  of  the  old  ones; 
and  the  new  axes  cannot  be  considered  physical  vectors,  because  their 
components  (i.e.  the  old  axes)  are  not  of  equal  physical  type.  Clearly, 
the  same  goes  for  every  point  in  such  an  axis  diagram,  except  for 
points  directly  on  one  of  the  axes.  Such  a  vector  could  be  considered 
of  "mathematical"  type,  since  it  may  not  represent  one  particular 
physical  quantity  at  one  particular  point  in  physical  space;  instead  it 
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may  represent  several  different  physical  quantities  (two,  if  the  axis 
diagram  is  two-dimensional)  not  necessarily  associated  with  one  and  the 
same  physical  point.  If  the  number  of  axes  (i.e.,  dimensionality)  on 
the  "diagram"  is  the  same  as  the  number  of  physical  quantities  involved 
in  the  system,  one  point  or  vector  in  the  diagram  will  describe  the 
state  of  the  whole  system,  i.e.,  such  a  vector  is  global. 

In  the  present  case,  the  mathematical  vector  is  represented  by 
a  physical  velocity  vector  field,  i.e.,  its  components  are  physical 
vectors  of  one  and  the  same  physical  type  (i.e.  velocity),  but 
associated  with  different  points  in  physical  space. 

A  small  random  perturbation  will  now  act  locally,  affecting 
a  certain  point  in  physical  space  and  its  neighbourhood  and  the 
physical  vectors  associated  therewith.  That  is,  such  a  perturbation 
can  affect  a  particular  subset  of  components  of  a  mathematical  vector, 
but  it  (almost)  always  will  affect  all  components  of  a  physical  vector 
at  once.  Such  a  perturbation  will  momentarily  uncouple  an  affected  point 
and  neighbourhood,  and  the  mathematical  vector  will  be  offset 
accordingly.  If  the  invariant  manifold  of  the  whole  system  is  a  3-torus, 
it  was  shown,  that  the  trajectory  of  this  mathematical  vector  can 
change  almost  randomly  under  these  perturbations.  Since  these  random 
perturbations  act  on  particular  points  or  areas  in  space  and  the 
trajectory  of  the  system  describes  its  development  in  time,  this  means 
nothing  else  than  uncoupling  in  space  leads  to  uncoupling  in  time. 

Clearly,  the  argument  for  a  3-torus  as  invariant  manifold  of  a 
system  to  cause  chaos  or  turbulence  can  be  derived  from  the  fact  that 
every  "physical"  space  is  three-dimensional  (i.e.,  3D  world).  In  an 
n-dimensional  physical  space,  the  projection  of  an  n-torus  would  be 
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necessary  to  create  the  type  of  structural  instability  against  random 

perturbations  that  would  lead  to  chaotic  behaviour.  At  the  present 

level  of  the  study,  virtually  nothing  specific  can  be  said  about  the 

requirements  on  random  perturbation  combinations,  which  can  cause  a 

drastic  change  in  the  state  space  trajectory.  However,  it  is  evident 

that  not  any  combination  of  random  perturbation  can  effect  such  a 

change  in  trajectory.  This  seems  to  relate  in  some  way  to  the 

requirements  Ruelle  and  Takens  impose  on  the  parameters  p^,p?,q^,q9, 

which  are  used  in  the  definition  of  their  vector  field  z  (see  page  17) 

and  the  fact  that  the  topology  specifying  the  type  of  neighbourhood 

between  the  vector  fields  z  and  to  cannot  be  any  "reasonable" 

2 

topology  but  must  be  C  . 

Another  interesting  question  in  connection  with  appropriate 
random  perturbation  combinations  would  be  how  the  probability  for  the 
occurrence  of  such  a  combination  changes  with  the  level  of  discretiza¬ 
tion  (coarser  of  finer)  and  the  type  of  probability  distributions  used 
on  the  random  number  generators . 

This  question  also  touches  the  problem  of  reality  of  the 
numerical  modelling  or  simulation  of  the  basic  theoretical  model.  As 
in  most  studies  involving  fluid  dynamic  numerical  simulation,  substantial 
compromises,  simplifications  and  approximations  are  forced  by  the 
complexity  of  the  problem  at  hand  (see  e.g.,  the  studies  of  the  cavity 
flow  problem  mentioned  in  Appendix  D  and  discretization  of  NS-equations 
in  Chapter  5) . 

One  of  the  biggest  simplifications  is  the  discretization  itself. 
Mathematically,  an  infinite-dimensional  function  space  is  approximated 
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by  a  finite-dimensional  one,  that  is,  by  a  set  of  discrete  functions. 

The  main  problem  here  is  a  drastic  change  of  compactness  properties, 
and  it  has  to  be  shown  for  a  particular  type  of  discretization,  that 
the  finite-dimensional  approximated  spaces  have  the  right  upper  limit, 
i.e.,  the  correct  infinite-dimensional  solution  space.  A  considerable 
amount  of  this  theoretical  groundwork  is  collected  in  Temam  [12]  .  On 
the  other  hand,  it  can  be  argued  on  physical  grounds  that  such  an 
approximation  of  an  infinite-dimensional  space  must  be  reasonable,  since 
in  a  fluid  domain,  each  point  with  its  associated  variables  such  as 
pressure,  temperature,  velocity,  etc.  can  be  considered  a  quite  accurate 
representative  for  a  certain  neighbourhood  around  it,  as  long  as  this 
neighbourhood  is  chosen  small  enough.  Therefore,  the  selection  of  a 
finite  number  of  points,  whose  neighbourhoods  define  a  partition  on 
the  fluid  domain  directly  relates  to  the  transition  from  an  infinite¬ 
dimensional  solution  space  to  a  more  or  less  accurate  finite-dimensional 
approximation,  depending  on  the  number  of  selected  points  (i.e. 
coarseness  or  level  of  discretization)  and,  therefore,  on  the  size  of 
their  neighbourhoods.  The  use  of  finite-element  techniques  with 
appropriate  overlapping  basis  functions  or  weight-functions  can  there¬ 
fore  be  considered  as  just  a  technical  trick  for  reducing  the  number 
of  selected  points  further,  without  losing  too  much  accuracy  while 
incorporating  certain  conditions,  such  as  incompressibility. 

A  further  concession  for  the  sake  of  technical  feasibility  has  to 
be  made  by  substituting  a  variational  formulation  for  the  original 
formulation  of  the  problem, as  a  consequence  of  the  employed  Galerkin 
method.  The  variational  formulation  might  permit  more  solutions  than 
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the  original  one. 

The  restriction  to  rather  specific,  i.e.,  Dirichlet  boundary 
conditions  on  a  bounded  subdomain  within  the  entire  fluid  domain  for 
the  time  dependent  perturbations  is  a  direct  consequence  of  the  rather 
high  generality  for  the  types  of  mean-flow  profiles.  These  profiles 
can  be  selected  without  any  considerations  for  their  underlying  boundary 
conditions.  The  motivation  for  this  approach  was  the  fact  that,  in 
experiment,  turbulent  "spots"  can  be  observed  in  an  otherwise  laminar 
environment,  although  in  the  present  model,  these  "spots"  have  to  be 
exactly  rectangular  due  to  the  type  of  space  discretization.  This 
drawback,  however,  could  be  reduced  by  using  finer  discretization  and 
identifying  the  rectangular  boundaries  for  the  time  dependent  perturbations 
with  the  boundaries  of  the  mean  flow  under  an  appropriate  flow  profile.* 
Another  important  factor,  which  requires  further  investigation, 
is,  as  already  mentioned,  the  level  of  discretization  and  the  shape 
and  size  of  the  underlying  rectangular  elements,  which  can  be  varied 
along  the  three  axes  of  the  domain  under  consideration.  This  would  vary 
the  level  of  discretization  within  the  domain  spatially,  and  thereby 
permit  to  focus  the  study  onto  the  more  critical  and  interesting  areas 
of  the  flow.  In  connection  with  this  local  and  global  varying  of  the 
discretization  level,  it  would  be  interesting  to  observe  how  the  critical 
value  of  the  bifurcation  parameter  and  other  output  parameters  of  the 
bifurcation  program  BIFOR 2  change,  especially  the  approximation  for 
and  solution  trajectory  for  . 

The  simple  geometric  configuration  for  such  a  study  would  be  the  one 
of  the  already  mentioned  cavity  flow  (see  Appendix  D) . 
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Also,  with  the  local  and  global  varying  of  the  discretizations, 
an  analysis  of  the  bifurcational  behaviour  of  different  mean  flow 
profiles  would  be  important.  Amongst  the  four  profiles  tested  in  this 
study  the  ones  with  the  more  interesting  results,  i.e..  Flows  III  and 
IV,  seem  to  be  too  complex  in  their  profiles  to  be  realistic.  That  is 
to  say,  simpler  profiles  might  show  a  more  interesting  bifurcational 
behaviour  for  a  higher  level  of  discretization.  The  main  deficiency  of 
the  numerical  simulation,  i.e.  stability  of  mean  flow  U  for  negative 
viscosities  only  (apart  from  the  case  of  exception  mentioned  in 
Appendix  C)  could  be  a  consequence  of  the  downgrading  induced  by  the 
discretization  as  described  previously,  such  as:  Dirichlet  boundary 
conditions,  variational  principle,  low  level  of  discretization.  However, 
the  main  reason  for  this  phenomena  is  most  likely  the  neglect  of  a 
pressure  perturbation  gradient.  That  would  mean  absence  of  a  driving 
force,  which  is  compensated  for  by  negative  viscosities  thereby 
creating  a  driving  force. 

However,  apart  from  this  negative  viscosity  and  despite  the  fact 
that,  for  reasons  of  economy,  the  lowest  possible  level  of  discretiza¬ 
tion,  i.e.,  4  points,  had  to  be  chosen,  the  results  show  in  general  a 
qualitatively  correct  behaviour  for  flow  IV.  The  trajectories  of 
V^  +  V^  clearly  show  an  increase  of  the  level  of  instability  of  V^ 
when  increasing  the  value  of  the  bifurcation  parameter  n  .  The 
physically  relevant  range  for  n  seems  to  be  .242  <_  n  <  .440,  i.e., 
between  the  critical  value  of  the  first  Hopf  bifurcation  and  the  value 
for  occurrence  of  unbounded  growth.  It  can  be  assumed  that  beyond 
n  a-  .449  the  shear  of  flow  IV  becomes  unrealistically  high  or  at 
least  too  high  for  the  low  level  of  discretization.  Moreover,  since 
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is  represented  by  its  first-order  approximation,  which  is  exact 
only  at  n  =  nc  =  .242,  its  inaccuracies  at  p  %  -449  might  be  too 
high  and  therefore  contribute  to  a  breakdown  of  the  model . 

The  most  interesting  and  most  promising  parameter-range,  as  far 
as  compliance  with  the  hypothesis  of  a  3-torus  as  invariant  manifold  is 
concerned,  seems  to  be  the  neighbourhood  of  p  .3.  The  unperturbed 
trajectory  projections  give  reason  to  suspect  the  presence  of  at  least 
two  frequencies,  which  are  close  to  each  other  and  non-rationally  related. 
Whether  this  is  only  a  transitional  stage  of  the  trajectories  between 
the  orbit  for  V  and  a  double-loop  limit  cycle,  as  can  be  seen  for 
trajectories  at  p  >  .33,  would  require  monitoring  over  considerably 
longer  time  intervals  as  could  be  done  with  the  CPU  time  available. 

The  fact  that  the  randomly  perturbed  trajectories  even  for  n>  -33  show 
no  tendency  to  approach  such  a  double- loop  limit-cycle  suggests  that 
if  such  a  limit- cycle  is  present  for  p  <  .33,  it  must  be  a  very  weak 
attractor  if  any  attractor  at  all. 

The  perturbed  spectra,  on  the  other  hand,  give  reason  to  assume 
the  range  of  turbulence  and,  therefore,  the  existence  of  a  3-torus 
as  invariant  manifold,  for  a  parameter  p  between  %  .33  and  at 
least  %  .445.  However,  this  has  to  be  looked  upon  with  a  certain 
amount  of  suspicion  and  doubt,  since,  as  already  mentioned  in  the 
discussion  of  the  spectra,  the  sample  sizes  had  to  be  kept  too  small 
to  permit  definite  conclusions  about  the  spectra  and  their  slopes.  For 
a  comparison  to  experimentally  obtained  spectra,  the  sample  size  would 

4 

have  to  be  increased  to  at  least  1.2 x  10  ,  which  is  about  the  lower 
limit  for  experimentally  obtained  samples.  Wind  measurements  for 
spectral  analysis  are  usually  done  over  periods  from  20  minutes  up  to 
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one  hour,  with  a  number  of  10  to  as  many  as  100  measurements  per  second 

and  subsample  sizes  covering  one  minute.  This  gives  overall  sample 

4  5 

sizes  from  1.2  x  10  up  to  3.6  x  10  .  For  smaller  sample  sizes 
smoothing  difficulties  are  encountered,  so  a  definite  spectral  slope 
cannot  be  determined.  Due  to  restrictions  on  available  CPU  time, 
the  numerical  analysis  in  general  and  the  spectral  analysis  in 
particular  of  this  study  could  not  be  advanced  beyond  its  present  level, 
which,  especially  in  comparison  with  the  orders  of  magnitude  of  sample 
sizes  for  experimentally  obtained  spectra,  cannot  be  considered 
anything  more  than  preliminary. 
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APPENDIX  A 


Some  of  the  concepts  used  in  Chapter  2  and  3  are  explained 
here.  Except  for  the  description  of  the  Cantor  set  and  the  strange 
attractor,  all  definitions  and  some  of  their  interpretations  are 
adopted  from  [9] : 

1.  THE  CANTOR  SET 
Construction 

The  best  known  way  of  constructing  a  Cantor 
set  is  to  trisect  the  closed  unit  interval  I  =  [0,1]  at  the  points 
1/3  and  2/3  and  then  delete  the  open  interval  (1/3,  2/3),  called 
the  'middle  third'.  C^  denotes  here  the  remainder  of  the  points 
in  I ,  i . e .  , 


C1  =  [0,1/3]  u  [2/3,  1] 

We  now  trisect  each  of  the  two  segments  in  C^  at  1/9  and  2/9  and 
7/9  and  8/9,  and  then  delete  the  'middle  third'  from  each  segment, 
i.e,  (1/9,  2/9)  and  (7/9,  8/9).  Let  C?  denote  the  remainder  of 
the  points  in  C^,  i.e., 

C2  =  [0,1/9]  u  [2/9, 1/3]  u  [2/3, 7/9]  u  [8/9,1] 


If  we  continue  in  this  manner  we  obtain  a  descending  sequence  of  sets 


C1  =  C2  =  c3  = 


where  C  consists  of  the  points  in  C  ,  excluding  the  'middle  thirds', 
m  m-1 


93 


94 


Observe  that  C  consists  of  2m  disjoint  closed  intervals  and,  if 
we  number  them  consecutively  from  left  to  right,  we  can  speak  of  the 
odd  or  even  intervals  in  C  .  The  Cantor  set  C  is  the  intersection 


of  these  sets,  i.e. 


m 


C  =  n{C.  :  i  e  N} 


Some  Properties 

1.  C  is  non  denumerable.  To  see  this,  define  a  function  f  on  C 
as  follows: 


f  (x)  =  <ax,a2, . . .  > 


where 


a 

m 


0  if  x  belongs  to  an  odd  interval  in 
2  if  x  belongs  to  an  even  interval  in  C 

m 


C  is  equivalent  to  the  set  of  sequences  <  a^ , a9 , .  .  .  >  ,  where 
a^  =  0  or  2,  which  has  cardinality  2a'*'e^^  nuH  which  is  the 
continuum.  Therefore  every  point  in  C  is  an  accumulation  point. 

2.  C  has  measure  zero.  The  measure  of  the  complement  of  C  relative 
to  I  =  [0,1],  i.e.,  the  union  of  the  middle  thirds,  equals 


12  4  8 

—  +  —  4*  -  +  -  + 

3  9  27  81 


...  =  1 


Since  the  measure  of  I  is  also  1,  the  measure  of  C  must  be  zero. 
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2.  THE  HORSESHOE  DIFFEOMORPHISM* 

The  present  introduction  is  rather  informal.  For  its  original 
introduction,  see  the  paper  by  Smale  [10],  who  discovered  this 
diffeomorphism  in  1966. 

2 

Take  a  rectangle  R  in  the  plane  ]R  .  Stretch  it  and  bend  it 
into  a  horseshoe  shape  as  shown  in  Figure  5  (vertices  A,B,C,D  go  to 
A'jB'jC'jD')  and  place  it  on  top  of  its  original  position.  This  can 

be  done  mathematically  using  smooth  but  possibly  nonanalytic  functions 

2  2 

pieced  together  to  obtain  the  proper  diffeomorphism  f  :  ]R  -*■  ]R  . 

It  is  important  that  this  diffeomorphism  stretches  the  two  rectangular 

strips  R  n  f  ^(R)  (consisting  of  points  of  R  which  remain  in  R 

after  f  is  applied)  linearly  by  a  factor  >1  in  the  vertical 

direction  (parallel  to  AD)  and  compresses  these  strips  by  a  factor  <  1 

in  the  horizontal  direction  (parallel  to  AB) . 

- 1  -2 

Now  look  at  R  n  f  (R)  n  f  (R) ,  i.e.,  the  set  of  points  which 

2 

remain  in  R  when  both  f  and  f  are  applied  to  them.  This  consists 
of  four  thin  rectangles,  two  in  each  rectangle  of  R  n  f  1 (R) . 

Similarly,  R  =  R  n  f  1 (R)  n  ...  n  f  n(R)  consists  of  2n  thin 
rectangles  contained  pairwise  in  the  2n  ^  rectangles  of  R  ^  ^ 

(see  Figure  6) . 

The  set 


R  =  R  n  f_1 (R)  n  f'2(R)  n  . . . , 

—  CO 

which  consists  of  those  points  of  R  which  remain  in  R  after  any 
number  of  positive  iterations  of  f,  is  the  cartesian  product  AB  x  C, 

k 

A  diffeomorphism  is  a  differentiable  mapping  with  differentiable  inverse. 
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Figure  5.  Horseshoe-Diffeomorphism. 

First  Mapping 


Figure  6.  Horseshoe-Diffeomorphism. 

Second  Mapping 

1 
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where  AB  is  the  closed  interval  corresponding  to  the  width  AB  of  R 
and  C  is  an  example  of  a  Cantor  set.  Clearly  the  diffeomorphism  f 
could  be  constructed  to  have  this  Cantor  set  develop  exactly  like  the 
one  constructed  by  taking  out  the  middle  thirds.  Likewise,  the  set 
consisting  of  those  points  of  R  which  remain  in  R  under  all  negative 
iterations  of  f  is  of  the  form  C  x  AB.  If  A  denotes  R  n  R 

—  CO  00 

then  we  see  that  A  is  a  closed  invariant  set  for  f,  meaning  that 
f  takes  A  onto  itself,  and  A  is  homeomorphic  to  C  x  C,  which  can 
be  shown  to  be  in  fact  homeomorphic  to  C.  Furthermore,  it  can  be 
shown  that: 

(1)  There  is  an  infinite  number  of  periodic  points  of  f  in  A  . 

(2)  There  are  points  which  have  orbits  that  are  dense  in  A  . 

(3)  The  effect  of  f  near  A  is  to  contract  in  a  horizontal  direction 
and  expand  in  a  vertical  direction. 

Clearly,  this  construction  can  be  modified  using  a  disk  instead 
of  rectangles  and  a  diffeomorphism  f  such  that  f(R)  c  R  (see  Figure 
7) .  This  is  the  variation  used  by  Ruelle  and  Takens  and  it  is  easy  to 
find  a  subset  homeomorphic  to  the  rectangle  ABCD,  which  develops  in 
the  way  described  above  (see  Figure  6) . 

3.  NONWANDERING  SETS 
Definition 

A  point  x  e  X  is  said  to  be  nonwandering  with  respect  to  a 
mapping  f  if  for  every  neighbourhood  U  of  x,  there  exists  an 
n  >  0  with  U  n  f_n(U)  ^  0  .  The  set  of  all  nonwandering  points  is 

called  the  nonwandeving  set.  That  means,  the  point  x  under  the 
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mapping  f  will  sooner  or  later  return  arbitrarily  close  to  any 

position  it  ever  takes  or  took  under  f  .  E.g.,  for  a  quasiperiodic 

k  k 

motion  on  the  k-torus  T  we  have  T  itself  as  the  nonwandering  set. 

4.  STRUCTURAL  STABILITY 

In  the  following  definitions  we  use  the  notations: 

p 

X  (M)  is  the  set  of  all  vector  fields  differentiable 

up  to  at  least  r^  order  defined  on  the  manifold  M. 

Diff  (M)  is  the  set  of  all  diffeomorphisms  differentiable 

th 

up  to  at  least  k  order  defined  on  the  manifold  M. 

Definition 

r  k 

Two  vector  fields  x,  y  e  X  (M)  are  C  equivalent  if  there 

k 

exists  h  e  Diff  (M)  mapping  x-traj ectories  to  y-traj ectories  (and 
vice  versa),  preserving  their  sense  of  motion,  but  not  necessarily 
their  parameterization  by  time. 

Definition 

A  vector  field  x  €  X  (M)  is  structurally  stable  if  there  is  a 

t  r 

neighbourhood  N  of  x  in  X  (M)  (N c x  (M))  such  that  every 
vector  field  in  N  is  C°  equivalent  (i.e.,  topologically  equivalent 
or  homeomorphic)  to  x.  This  means  that  the  (topological)  structure 
of  the  trajectories  of  x  does  not  change  if  x  is  slightly  perturbed 
so  that  it  becomes  some  other  vector  field  x',  say,  which  is  still  an 
element  of  the  neighbourhood  N.  In  the  most  simple  terms:  if  any 
vector  field  close  enough  to  x  has  trajectories  very  much  like  those 
of  x,  then  the  vector  field  x  is  structurally  stable. 
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5.  MORSE-SMALE  SYSTEMS 
Definition 

A  fixed  point  p  of  a  vector  field  x  is  called  hyp er>b olio ,  if 
the  linearized  vector  field  at  p  has  no  purely  imaginary  or  zero 
eigenvalue.  That  means  the  eigenvalues  of  the  linearized  field  x 
at  p  must  have  either  positive  or  negative  real  parts, respectively . 
Therefore,  the  vector  field  x  in  a  certain  neighbourhood  U  of  p 
must  have  components  either  radially  directed  outwards  or  inwardSj 
respectively,  w.r.t.  p. 

The  set  of  all  trajectories  of  a  vector  field  is  sometimes  called 
the  flow  and  denoted  by  <j>  .  More  specifically,  cf>  (p)  is  used  to  denote 
the  point  up  to  which  the  point  p  "travels"  within  the  time  t  under 
the  influence  of  the  vector  field.  Looking  at  the  obvious  properties 

<f>t  *  ‘fr-t  =  ^o  =  identity 
ot  •  *s)  •  *u  =  *t  •  0f>s  •  *u)  =  *t+s+u  > 

it  is  evident  that  cj>  is  a  commutative  Lie  group. 

Definition 

Let  <j)  be  the  flow  of  a  vector  field  x  on  a  manifold  M  and 
p  a  hyperbolic  fixed  point.  The  set 

W  ±  (p)  =  {q  e  M  |  <f>  (q)  ->  p  for  t  ->  ±°o} 

is  called  the  stable  (for  +)  or  unstable  (for  -)  manifold  for  p. 

Clearly,  a  point  on  the  stable  manifold  for  p  is  attracted 
towards  p  through  the  vector  field  x  for  positive  times,  whereas 
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a  point  on  the  unstable  one  is  attracted  to  p  for  negative  times. 

Definition 

r 

A  vector  field  x  e  X  (M)  ,  is  Movse-Smale  (MS)  if: 

(1)  The  number  of  fixed  points  and  periodic  orbits  is  finite,  and 
each  is  hyperbolic. 

(2)  The  set  of  nonwandering  points  consists  of  fixed  points  and  periodic 
orbits  only. 

(3)  All  stable  and  unstable  manifolds  intersect  transversally . 

For  two-dimensional  manifolds  the  following  statements  hold: 

1.  structural  stability  implies  MS. 

2.  MS  systems  are  dense  in  Diff^  (M)  or  X'*'  (M)  . 

3.  Structurally  stable  systems  are  dense  in  Diff  ^(M)  or  X  ^(M) . 

For  manifolds  with  dimensions  higher  than  2  only  MS  =>  structural 
stability  has  been  proved  by  Palis  and  Smale.  The  rest  of  the 
statements  is  false  in  general. 

6.  GENERICITY 
Definition 

Let  S  be  any  topological  space.  A  subset  G  of  S  is  called 
vesi-dual  if  it  is  the  intersection  of  a  countable  number  of  sets,  each 
of  which  is  both  open  and  dense  in  S. 

The  motivation  behind  this  definition  is  to  distinguish  between 
two  different  types  of  denseness.  For  example,  both  the  set  of  rational 
numbers  Q  and  the  set  of  irrational  numbers  R  are  dense  in  the  set 
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of  real  numbers  ]R  .  However,  n  is  in  the  same  sense  "denser"  in 
]R  then  Q;  there  are  more  irrationals  than  rationals.  Since  the 
rationals  are  countable: 


Q  =  (qrq2 


one  can  express  the  irrationals  as  the  countable  intersection: 


cx> 


(*-  {qn})- 


n  =  n 

n=l 


This  construction  of  the  irrationals  follows  the  definition;  they  are 
therefore  a  residual  subset  of  the  reals,  whereas  the  rationals  are  not. 

The  definition  of  genericity  is  just  a  slight  generalization  of 
the  previous  one: 

Definition 


A  property  p  of  elements  e  of  a  set  S*  is  called  genev-ia  if 


these  elements  form  a  set  P  which  contains  a  residual  subset  of  S. 

Replacing  a  real  number  by  another  real  by  adding  an  arbitrary 
"perturbation",  the  probability  of  the  new  real  to  be  rational  is  zero 
whereas  the  probability  to  be  irrational  is  one.  That  is,  the  generic 
property  is  adopted.  This  shows  how  the  concept  of  genericity  can  be 
understood  and  applied  in  connection  with  the  vector  fields  discussed 
at  the  end  of  Chapter  2. 


In  general,  S  will  contain  elements  other  than  the  e  as  well. 


APPENDIX  B 


I.  THE  DIRAC  A-FUNCTION 

The  A-function  can  be.  defined  as  any  distribution  with  the 
following  properties: 

(i)  A  (x)  =  0  for  x  ^  r 

(ii)  /  Ar(x)dx  =  1 

(iii)  /  Ar(x)f(x)dx  =  f  (r)  . 


The  integration  can  go  over  more  than  one  dimension  and  the  range  of 
integration  naturally  has  to  include  the  point  r. 

Analytic  representations  of  A^(x)  are  e.g.: 


ca 

A  (x)  =  (2tt)  /  exp  [ico  (x-r)  ]  dm 

—  00 


or: 


A 

r 


(x)  =  lim 


sin (fix) 

TTX 


In  the  present  context,  the  property  (iii)  is  used  with  the 
following  notation: 


<Ar(x)  I  f(x)>=  f(r), 

i.e.,  it  is  written  as  a  scalar  product  of  two  vectors  defined  in 
function  space. 
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2.  CONSTRUCTION  OF  A  .LINEAR  10  x  10  SYSTEM  FOR  THE  DETERMINATION  OF 
THE  COEFFICIENTS  IN  THE  BASIS  FUNCTIONS 

As  noted  in  the  text,  the  basis  functions  W^(x),  r  e  S  are 
polynomials : 


W 


(x)  =  c  +  y  c.x.  +  y  c. .x.x. 

v  11  . u  .  Ill] 

l  1<J  J  J 


CB1) 


x  =  (x1,x2,x3) 


With  the  following  conditions: 


Four  Nodal  Value  Conditions 


W  (r)  =  1 
r 


(B2) 


Wr(rL)=  0,  L  =  1,2,3 


Here  the  r  denote  the  three  centroids  on  the  other  three  faces  of 

1j 

the  tetrahedron. 


Three  Conditions  for  Zero  Divergence  Requirement 


8W  (X) 

dXldX2dX3  *  °-  1  *  1'2’3 


(B3) 


(integration  over  one  tetrahedron) 


Three  Gradient  Conditions 


3W  (X) 
r 

ax. 

i 


=  0 ,  i  =  1 , 2 , 3 


(B4) 


X=r 
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To  apply  these  conditions,  the  basis  functions  W(X)*  have  to 
be  transformed  into  barycentric  coordinates. 

Notations : 


X  with  indices  i,j,k  denotes  cartesian  coordinates 
with  indices  £,m,n,o  denotes  barycentric  coordinates 
X  with  double  indices  i,j,k;  £,m,n,o  denotes  cartesian  coordinates 

(i,j,k;  first  index)  of  the  four  tetrahedral  vertices  (£,m,n,o;  second 
index) 

For  the  relations  between  cartesian  and  barycentric  coordinates, 
we  write 


or : 


X11 

X21 

X31 

X41 

X12 

X22 

X32 

X42 

X13 

X23 

X33 

X43 

1 

1 

1 

1 

X. 

1 


X.  rxT 
lL  L 


(B5) 


(B6) 


Using  (B6) ,  the  polynomial  (Bl)  and  its  derivatives 


aw 

8Xk 


(X) 


=  c. 


C..  X. 
lk  1 


C,  ,  X.  , 
kk  k 


k  =  1,2,3 


Subscript  r  dropped  for  conciseness. 


(B7) 
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can  now  be  transformed  directly: 


3 

"4 

"4 

'4 

W[X(A)]  =  C  +  I  C 

1 

i  ,“>‘J 

+  l  c.  . 

•  •  1] 

j 

- 1 

j 

x 

•r— i 
•H 

X 

[X-J  i — ! 
_ 1 

CB8) 


3W 

3X, 


[X(A)]  = 


3 

I 


Ck  +  4-  "ik 
1 


'4 

"  4 

l  XiLXL 

_  Li 

+  ckk 

[  hih 

(B9) 


Relations  (B5)  or  (B6) ,  respectively  show  linear  dependence 
between  the  barycentric  coordinates  A  .  Therefore,  in  order  to 

lj 

perform  the  integrations  for  the  zero-divergence  condition,  one  of  the 


A^  must  be  eliminated.  The  vertex  associated  with  the  eliminated 

barycentric  coordinate  is  the  origin  of  a  skew  coordinate  system 

spanned  by  the  three  edges  of  this  vertex. 

Taking  the  vertex  X  ,  i  =  1,2,3  as  origin,  the  integral  of 

1  . 

9  W 

—j—  over  the  tetrahedron  may  be  written  as 

o  A, 

k 


1  ^t  8Xk  tX(-Xl,X2,X3^ 


3  X . 

_ i_ 

3  At 


dbdX2dX3 


(BIO) 


Here 


3Xt 


is  the  functional  determinant  with  L  f  4 


Eliminating  A^  from  (B6)  gives: 


Xi  ■  I  hl\  +  W1' VW 

Lj 


=  V  (x.T  -  x.,)at  +  x., 

£  lL  i4  L  i4 


(Bll) 


therefore : 
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3X. 

X11 

-X14 

X21 

-X24 

X31 

-X34 

i 

_ 

x 

-  X 

X 

-  X. 

X  „ 

-  X 

3Al 

12 

14 

22 

24 

32 

34 

L  LM  - 

X,  „ 

-  x,  „ 

x„. 

-  x_ 

x_ 

-  X 

13 

14 

23 

24 

33 

34 

[ X . t  -  X.  ,] 
lL  i4  dgt 


(B9),  (Bll),  (B12)  into 

(BIO): 

3 

'3 

I  =  [X..-X.J,  / 

lL  i4Jdet 

V?  Cik 

l 

l  (XiL-Xi4^L  +  Xi4 

__  Li 

+  C 


kk 


T  (X.T  -X.  JX_  -  X.  . 
“  lL  i4  L  i4 

L 


dAldA2dA3 


3 

"3 

L  lL  i4Jdet 

1  c,  + 1  C., 

o  k  h  lk 

l 

y  (X..-X.  JIT  +  I  X.  . 

“  lL  i4^  L  o  i4 

+  C 


kk 


T  (X..-X.JI.  +  I  X.  . 
f  ^  lL  i4J  L  o  i4 


with: 


I  =  /  dX^dX-dX.  =  i 
o  J  _  12  3  6 


1 L  ALdAldA2dA3  24 


The  integrals  Iq,I^  are  most  easily  evaluated  by  use  of 


Hammer's  formula: 


pl!p2!p3! 


/  X^X^X^dX  dX  dX_  =  ■■■■  --J)- - 

JT  1  2  3  122  (p1+p2+p3+3) 


CB12) 


(B13) 


CB14) 


(B15) 


The  equations  for  the  ten  conditions 


(B2) ,  (B3) ,  (B4)  can  now  be 
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written  down  explicitly: 


4  Nodal  Value  Conditions: 


W(r4)  =  1  : 


is  the  centroid  on  the  face  opposite  to  the  barycentric 

origin  (x14 >X24’X34) •  Using  (B8)  and  r4  (A^A^A^  =  (1/3 , 1/3 , 1/3)  ; 

A  =  0  : 

4 


1  = 


3 

“  3 

3 

‘  3 

",  3 

=  c  +  y  c. 

t  y  x .  T 

+  y  c. . 

\  y  x. _ 

T  J  X  . 

b  l 
l 

[3  L  iLJ 

.  ij 

i<j 

_3  L  lL^ 

L3  L  jLJ 

(B16) 


W(r  )  =  0  : 


Coordinates  of  r^  :  A^  =  A2  =  A4  =  y  ,  A^  =  0 


0  =  C  +  l  C: 


—  y  x. 

3 

+  y 

C.  . 

—  y  x 

4  y  x.T 

L3  L=1 , 2,4  lLJ 

L 

i<j 

ij 

3  L=1 , 2 , 4  lL 

L=l,2,4  jL_ 

(B17) 


W(r  )  =  0  : 


Coordinates  of  r2  :  A^  =  A^  =  A4  =  y  ,  A^  =  0 


o  =  c  +  y  C. 

V  1 
1 


W(r  )  =  0  : 


—  y  x. 

3 

+  y  c. . 

i^ 

■—i 

X 

_ 1 

t  y  x .  t 

L3  L=l,3,4  lLJ 

•  •  ij 

i<j  J 

_  L=1 ,3,4  lL_ 

3  L=l,3,4  ^ 

(B18) 


Coordinates  of  r^  *  ^2  =  X3  =  ^4  =  T  ’  Xi  =  0 
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0  = 


=  C  +  l  C. 
4  1 


l  y  x 

3 

.  y  r 

1  y  x 

1  y  x 

3  l=2,354  iL 

ill 

3  L=2 ,3,4  iL 

3  L=2,3,4  ^'L 

(B19) 


3  Zero  Gradient  Conditions: 

3W(r  ) 

Using  Equation  (B9)  (k  =  1,2,3)  and  — — -  =  0 

<3  A, 


0  =  j-  c. 

6  k 


3 

'i  3 

"l  5 

+  y  c.. 

—  y  x 

+  C  , 

—  y  x. 

4  ik 
i 

L3  L  iLJ 

kk 

L3  L  iLJ 

k  =  1,2,3  (B20) 


3  Zero  Divergence  Conditions  over  Tetrahedron: 

Using  Equation  (B13)  (k  =  1,2,3)  and  1=0: 


1  3 

0=^0.  +  {  C.,  ... 
6  k  4  lk  \  24 

l 


>4  £  XiL^  +  Ckk^24  £  XiL^ 


k  =  1,2,3  (B21) 


introducing  the  technical  substitution: 


S. 

l 


4 

l  XiL 

L  1L 


the  linear  system  for  the  coefficients  can  be  written  down, 


1  3  1  ^ 

1)  0  =  C  +  4-  y  C.  [S.-X.J  +  i  y  C.  .  [S.-X.J  [S.-X.J 
3  4  iL  i  ilJ  9  .L  .  ij  L  l  ilJ  L  j  jl 


1<J 

3 


2)  0  =  C  +  ~  y  C.  [S.-X.J  +  4  y  c.  .  [S.  -X.  „]  [S.-X.J 
3  4  iL  i  i2  9  .L .  ij  i  i2J  j  j2J 


i<J 

3 


3)  0  -  c  ♦  J  ci[si-xi3]  J  I  C  [S.-Xi3][S  -X  ] 

1  1<J 

3  3 

4)  1  =  C  +  4  y  C.  [S.-X.J  +  4  y  C.  .  [S.-X.J  [S.-X.J  . 
3  4  iL  l  i4J  9  .LJ.  ij  l  i4  j  j 4  y 


i<] 


\ 


Nodal  Values 
Conditions 
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5)  0 

6)  0 

7)  0 

8)  0 

9)  0 

10)  0 


C1  + 


C2  + 


C3  + 


—  C 
6  L1 

—  C 
6  L2 


—  C 
6  3 


I  ClltSl-Xl4J  +  T  WV 

+  J  C12[SrX14] 
J  C33tS3'X34^  +  J  C23^S2'X24^ 


+  jSDW 

+  J  C23^S3'X34-I 

+  i  c3i[srxi4] 


Zero  Gradient 
Condition 


+  12  C11S1  +  24  C 1 2 S 2  +  24"  C31S3 
+  12  C22S2  +  24  C 1 2 S 1  +  24  C23S3 
+  12  C33S3  +  24  C23S2  +  24  C31S1 


Zero  Divergence 
Condition 


(B22) 

The  system  is  written  down  for  the  basis  function  W  (X),  i.e., 

r4 

the  one  which  takes  value  1  at  the  centroid  in  the  face  opposite  to 

vertex  No.  4.  It  can  be  converted  to  the  systems  for  W  (X),  W  (X), 

rl  r2 

W  (X)  by  replacing  the  seond  index  (i.e.,  4)  on  X  in  the  zero 
r3 

gradient  conditions  by  1,2,3  respectively  and  by  permuting  the  1  on 
the  left  hand  side  of  Equation  (4)  in  the  system  to  the  left  hand  side 
of  Equation  1) ,  2) ,  3) ,  respectively. 


3.  THE  PROGRAM  DSCRB 

DSCRB  discretizes  a  given  three-dimensional  rectangular  domain  D 
into  tetrahedrons,  calculates  the  ten  coefficients  for  the  basis 
functions  of  all  nodes  according  to  the  linear  system  derived  in 
Section  2  and  uses  them  in  the  basis  functions  to  calculate  the  matrix 
elements  for  the  interior  nodes*  in  Equation  (5.14)  for  a  given  flow  U 
of  the  structure  (7.3).  Finally,  it  inverts  the  matrix  on  the  left 

•k 

To  induce  Dir ichlet  boundary  conditions,  the  nodes  on  3D,  consisting  of 
the  surface  planes  of  D,  are  omitted. 
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hand  side  over  to  the  right  hand  side  to  obtain  Equation  (5.20). 

The  program  ADV  complements  DSCRB  by  reassembling  the 
nonlinear  matrix  elements  WV  calculated  by  DSCRB  into  the 
additional  advective  matrix  WV  +  WV  in  the  Equations  (5.18), 

(5.19)  for  the  second  and  third  perturbations  and  and 

multiplying  it  with  the  inverse  of  the  left  hand  side  matrix  W 

on 

3 . 1  Discretization  of  the  Domain  D  in  the  Main  Line  Program* 

The  program  begins  by  partitioning  a  rectangular  domain  D  into 
blocks.  The  size  of  these  blocks  and  D  is  naturally  induced  by  three 
finite  sets  of  user-supplied  discrete  values  along  the  three  coordinate 
axes  (XIN,  YIN,  ZIN) .  At  the  same  time,  each  one  of  these  blocks  is 
split  into  five  tetrahedrons  according  to  configuration  A  or  B  (see 
Figures  12a, b) .  These  two  configurations  are  alternating  between  two 
blocks  sharing  a  face  in  order  to  obtain  an  admissible  discretization 
(see  Temam  [12],  p.  73),  which  means  that  directly  neighbouring 
tetrahedrons  have  to  have  common  faces  and  edges.  The  splitting  of  the 
individual  blocks  is  simply  done  by  reassigning  the  block  corners  as 
corners  to  the  five  tetrahedrons,  which  are  numbered  at  the  same  time 
from  one  to  five.  Next,  for  each  individual  tetrahedron  the  nodes  (i.e., 
centroids  of  the  faces)  are  identified  by  their  opposite  vertices,  the 
number  of  the  tetrahedron  and  the  coordinates  of  the  block  (identified 
by  the  three  numbers  of  discrete  values  along  the  three  axes,  which 
determine  the  left  bottom-front  corner  of  the  block) .  To  determine  row 
and  column  numbers  for  the  linear  matrices,  the  nodes  are  being  counted 

through  as  they  are  being  identified.  At  this  stage  each  node  is 

* 

See  also  3.10  for  details. 


* 
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checked  using  its  identification  whether  it  has  already  been  counted  in 
a  previously  encountered  tetrahedron  or  whether  it  is  located  in  3D. 

If  this  is  the  case,  the  node  is  skipped  by  the  count.  This  is 
possible  since  every  node  is  shared  by  two  tetrahedrons  (except  nodes 
in  3D)  and  carries  therefore  two  identifications.  Once  the  program 
established  one  identification,  it  can  be  used  to  determine  the  other 
one  and  thereby  whether  it  was  previously  encountered.  After  this 
counting  and  checking  procedure,  the  linear  matrix  elements  are  calculated 
by  subroutines  and  indexed  according  to  the  count  numbers  of  the  nodes 
and  their  place  in  the  fluid  space  subdecompositions. 

3 . 2  Matrix  Elements  for  Advective  Terms 

Before  the  matrix  generation  of  the  program  can  be  described,  it 
is  necessary  to  analyse  the  particular  structure  of  the  advective  matrix 


elements  W^7 
on 

and 

wV:U 

on 

(see 

(5.14)  and  (5 

.15)) 

for  a  mean  flow  U 

of  the  type 

(7.3). 

First , 

KW 

on 

and 

WV:U 

on 

can  be  combined 

into 

one  matrix  element: 

wuv  ♦ 

on 

WV:U 

on 

=  <Wr  (X)I  |Ad> 
0 

(B23) 

with 

Aa 

=  (U[V  :  Wr 

(X)]  +  [V  :  U]Wr 

(X)) 

(B24) 

as  the  discretized  version  of  the  advective  operator  defined  in  (4.3). 
Written  out  in  components  by  use  of  (5.10a),  (5.10b)  A  looks  like:* 


k 

W  (X)  is  simply  written  as  W  for  reasons  of  conciseness. 
rn 
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Ad  = 


W 


W 


w 


?U1  -  -*• 
gjT  +  [U  • 

w 

3ui 

W2 

w 

9U1 

3X3 

au0 

au 

au 

axx 

w 

ax7  +  [u  ■  vw] 

w 

au 
_ 2. 

w 

3U3 

w 

3U3 

ax 

3X2 

3X3 

+  [u  •  vw] 


(B . 25) 


with: 


U  •  VW  =  U,  +  U„  4r-  +  U  8W 


1  3X, 


2  ax2  '  ~3  ax3 


(B . 26) 


Substituting  the  structure  of  (7.3)  for  U  gives: 


Ad  -  ClMl  ♦  C2M2  ♦  C3M.  ♦  CUMU  ♦  C22M22  ♦ 

C12M12  +  C21M21  +  C23M23  +  C32M32  +  °1 2M13  +  C31M31 


with : 


(B . 27) 


M,  = 


8W 

3X, 


1 

h-> 

O 

O 

_ 1 

aw 

"1  0  0“ 

aw 

i 

o 

o 

_ 1 

0  1  0 

'  M2  =  3X 

0  1  0 

•  M3  =  3X 

0  1  0 

Lo  0  1 J 

Lo  0  lj 

3 

L°  °  lj 

M 


11 


F  -1  E 
11  11  3W 

WEnxl  +  xi  sT 


1 


li  aw 


i  ax 


o 


l 


o 


li  aw 


i  ax 


l 


o 
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M 


22 


X 


'22  aw 


2  ax. 


0 


0 


E22~1  E22  3W 

WE22X2  X2  3X„ 


X 


0 


'22  aw 


2  ax. 


M 


33 


'33  aw 


3  ax. 


0 


0 


X 


'33  aw 


3  ax. 


0 


E33_1  E33  aw 

WE33X3  +  X3  axT 


M 


12 


'12  aw 


2  ax. 


E12"1 

E12X2  *  W 


'12 


X 


aw 

ax. 


x 


'12  aw 


2  ax. 


M 


21 


'21  aw 


l  ax 

E21_1 

E21X1  '  W 


0 


0 


'21  aw 


Li  ax. 


o 


'2i  aw 


Li  ax. 
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M 


23 


X 


'23  3W 


3  3X, 


0 


0 


X 


'23  3  W 


3X 


0 


0 


E23-1 

E23X3  *W 


'23  3W 


3  3X, 


M 


32 


X 


'32  3W 


2  3X. 


0 


X 


0 


'32  3W 


2  3X. 


E32X 


E32_1 


W 


0 


'32  3  W 


2  3X. 


M 


13 


X 


'13  3W 


3X. 


0 


X 


13  3W 


3  3X, 


E13-1 

E13X3  *W 


0 


X 


13  3W 


3  3X. 


M  = 
31 


'31  3W 


3X. 


0 


E3l~1 

E31X1  *  W 


'31  3W 


1  3X. 


0 


'31  3W 


1  3X, 


This  rather  simple  structure  of  is  one  of  the  reasons  for 

the  particular  choice  of  a  flow  U  of  the  type  (7.3). 


j. 
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3 . 3  Subroutine  SORT 2 

For  each  individual  tetrahedron,  each  one  of  the  four  linear 
10  x  10  systems  for  the  coefficients  of  the  four  basis  functions  is 
assembled  according  to  (B22)  and  then  solved  by  the  function  sub¬ 
routine  MATINV.*  After  that,  the  coefficients  of  their  first  and 
second  partial  derivatives  are  determined. 

3 . 4  Function  Subroutine  QMLTI1 

Because  of  the  structure  of  U  (see  7.3)  and  the  basis  functions 
W  (X) ,  the  matrix  elements  are  linear  combinations  of  integrals  like: 


,  Pi  P2  P3 

J  X1  X2  VdXldX2dX3 


=  [X.T -X.  .]  ,  _  / 
L  lL  i4Jdet  J 


('X1L"X14')XL+X14 


l  (X2l~X24'IXL+X24 


^X3L  X34)XL+X34 


dX1dA2dA3 


(B28) 


with  integers  p^,  p^,  p^  . 

The  right  hand  side  is  the  transformation  into  barycentric 
coordinates  (using  (Bll) ,  (B12)).  The  determinant  ^XiL_Xi4^det 
evaluated  in  the  main  line  and  QMLTI1  evaluates  the  integral 
analytically  for  p^  +  p2  +  p^  <_  6**  by  use  of  Hammer's  formula  (B15)  . 


Supplied  by  its  author  D.  Oracheski,  Atmospheric  Environment  Service, 
Edmonton. 

•k  k 

A  version  for  Pi  +  P2  +  P3  —  ^  exists  as  well.  It  permits 

spatial  dependencies  in  U  upto  power  6.  However,  this  leads  to 
prohibitively  high  CPU  times. 
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This  is  the  main  reason  for  the  particular  structure  chosen  for  U  : 
most  spatial  dependencies  other  than  simple  power  laws  would  require 
numerical  integrations  or  might  have  very  little  variability  after 
discretization  and  integration.  This  would  restrict  the  choice  of 
possible  bifurcation  parameters. 

QMLTI1  consists  of  six  nested  loops  with  four  iterations  each 
for  the  four  summands  (X.T  -X..)AT  ,i  =  1,2,3,  X..  in  each  one  of 

the  up  to  six  factors  in  the  integral  of  the  right  hand  side  of  (B28) . 
That  means,  the  number  of  factors  (i.e.,  p^  +  p0  +  p^)  is  equal  to 
the  number  of  loops  used  and  the  remaining  inner  loops  are  skipped. 

The  computation  of  the  integrals  over  powers  of  the  barycentric 
coordinates  A.^,A2,A3  (i.e,  evaluation  by  Hammer's  formula)  is  done 

in  the  innermost  used  loop  only,  and  the  outer  loops  are  used  only  to 
assemble  the  correct  linear  combinations  of  these  integrals  with  the 
(X.T  -X..),  L  =  1,2,3,  X..  as  coefficients. 

Whenever  the  main  line  enters  a  new  tetrahedron,  QMLTIl 
evaluates  all  possible  integrals  satisfying  p^  +  p2  +  p^  6  on 
the  right  hand  side  of  (B28) .  After  the  results  have  been  multiplied 

with  t^iL-Xi4^det  ’  t^ieX  are  stored  in  the  array  QT  and  indexed 
according  to  their  exponent  combinations.*  This  method  is  considerably 
faster  than  recalculating  these  integrals  whenever  they  are  needed. 

3 . 5  Function  Subroutine  QMLTI2 

This  routine  uses  the  output  of  QMLTIl  and  S0RT2  to 

* 

This  storage  method  requires  an  offset  of  one,  i.e.: 

.  Px  P2  P3 

/T  X1  X2  X3  dX1dX2dX3  -  QT  (p1+l ,p?+l ,p3+l) ,  since  a  FORTRAN-66  array 
cannot  be  indexed  by  zero. 


•* 
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assemble  the  actual  matrix  elements.  It  consists  of  three  nested 

loops  with  variable  numbers  of  iterations  (and  variable  initial 

value  on  the  outer  loop) .  This  way,  multiplications  of  up  to  three 

factors  (i.e.,  nonlinear  matrix  elements  w'7  )  consisting  of 

onm 

polynomials  with  varying  numbers  of  terms  can  be  performed. 

The  proper  selection  of  exponent  combinations  p^p^jP^  from 
the  array  QT  is  controlled  by  the  data  sets  IEX,  IEY,  IEZ.  The 
basis  functions  W^(X)  are  used  throughout  in  the  following  order 
of  terms: 

Wr(X)  =  C(l)  +  C(2)Xj  +  C(3)X2  +  C(4)X3  +  C(5)X^  +  C(6)X2  +  C(7)X‘ 

+  C(8)X1X2  +  C(9)X7X3  +  CflOJXjX,  (B29) 

i.e.,  first  coefficient:  all  exponents  zero 

second  coefficient:  X^-exponent  =  1,  X? ,X^-exponents  =  0 


tenth  coefficient:  X^ -exponent  =  X?-exponent  =  1, 

X^- exponent  =  0. 

This  sequence  of  exponent  combinations  is  kept  for  the  first  and 
second  partial  derivatives  of  W  (X)  as  well,  where  it  goes  up  to 
fourth  or  first  coefficient  respectively  only.  IEX,  IEY,  IEZ  now 
contain  the  partial  sequences  for  the  X^,X  ,X^  exponents  which  are 
in  their  combination  the  above  sequence.  The  integers  IPX,  IPY,  IPZ 
are  either  used  to  provide  the  required  offset  on  QT  (see  last  foot¬ 
note)  or  to  increase  exponents  according  to  spatial  dependencies  of  U 
of  powers  higher  than  one. 
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3.6  Subroutine  QMTRXB 

QMTRXB  assembles  all  linear  3  x  3  matrices  on  fluid  space 

level  with  the  structure  (scalar)  x  (unit  matrix).  These  are: 

W  (LHS, (5. 13)) ,  WA  (viscosity  term,  (5.16)),  <W  I M..  >  ,  <  W  I M„ >  , 

o  o 

<W  |M^>  (advective  terms,  (B28)).  The  subroutine  works  straight 
o 

forwardly  using  QMLTI2  to  calculate  the  scalar  factor  of  these 
matrices  and  then  determining  the  row  and  column  numbers  of  these 
elements  in  the  final  3N  x  3N  matrices  by  the  formula: 

N  =  3(1  -  1)  +  J  (B30) 

with:  N  =  row  or  column  number  in  3N  x  3N  matrix 

I  =  count  number  of  node  rQ  (for  rows)  or  r^  (for  columns) 
J  =  row  or  column  number  in  3x3  matrices  on  fluid  space 
level . 

The  elements  are  indexed  by  these  row  and  column  numbers  and 

stored  accordingly  in  3N  x  3N  arrays.  Since  all  diagonal  matrix 

elements  such  as  e.g.,  W  consist  of  a  sum  of  two  integrals  over 

two  tetrahedrons  showing  a  face  with  r^  as  centroid,  whereas  all 

nondiagonal  elements  consist  of  integrals  over  one  tetrahedron  only 

(as  is  evident  from  the  fact  that  the  support  area  for  the  basis 

functions  W  (X)  consists  of  two  tetrahedrons  as  described  above) , 
n 

the  diagonal  elements  are  calculated  by  the  program  in  two  different 
parts.  To  have  these  two  parts  properly  added,  the  3N  x  3N  arrays 
are  initialized  to  zero  at  the  beginning  and  the  matrix  elements  are 


added  into  them. 
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3 . 7  Subroutine  QMTRXA 

QMTRXA  assembles  and  stores  all  linear  3x3  matrices  on 
fluid  space  level  with  a  structure  different  than  the  one  of  elements 
assembled  by  QMTRXB.  These  are  all  elements  involving  the  advective 
matrices  ^  listed  in  (B27) .  They  can  be  written  in  a 

general  form: 


M.  . 
ij 


X.30  [I,]  +  E.  .X.1*  •  W[I.  .] 

J  3X.  3J  ij  j  L  ijJ 


(B31) 


with  [I.  .1  as  a  3x3  matrix  with  I  =0  for  n  M;  m  /  i  and 
L  ijJ  nm 

I.  .  =  1. 
ij 

The  routine  consists  of  two  nested  loops  for  the  two  indices  i 
(outer  loop:  components  of  U  )  and  j  (inner  loop:  spatial  dependencies 
of  U  )  .  Again  QMLTI2  is  used  to  calculate  the  two  matrix  elements 
relating  to  the  two  different  terms  in  (B31)  separately  and  the 
exponents  are  user-supplied  determining  the  spatial  dependencies 

of  the  main  flow  U  as  specified  in  Equation  (7.3).  A  formula  of 
the  type  (B30)  calculates  rows  and  columns  for  the  nondiagonal  matrix 
elements  (i  j)  of  the  second  term  in  (B31)  and  they  are  added 
into  arrays  in  the  manner  described  in  Section  3.6. 

A  third  loop  inside  the  two  i,j -loops  distributes  all 
diagonal  elements  (i=j)  into  storage  arrays,  according  to  a  formula 
like  (B30)  . 

3.8  Subroutine  QMTRXN 

As  is  evident  from  the  3x9  subdecomposition  of  the  nonlinear 
matrix  elements  in  (5.12a),  the  nine-dimensional  product  vector 


N< 
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V, . (r  )V, . (r  )  can  be  contracted  into  a  six-dimensional  product  vector 
lj  n  lj  m  ^ 

for  n  =  m,  since  equal  components  of  V(r  )  become  commuting  and 
therefore: 


al<Vli<VVli<rn»  +  a2<VUCrn)Vli(rn»  ’  tal  +  a2^VliCVVli(rn”’  (B32) 


and  the  3x9  subdecomposition  can  be  contracted  into  a  3x6 
subdecomposition  as  well: 


3W 

3X1 

0 

0 


0 

3W 

3X2 

0 


0 

0 

3W 

3X_ 


3W 

3X2 

3W 

3X, 


0 


3W 

9X3 

0 

3W 

3X, 


0 

3W 

3X3 

3W 

3X~ 


(B33) 


3  W 


m 


QMTRXN  uses  QMLTI  2  to  calculate  the  three  elements  <  |  |- 

o  n  i 

i  =  1,2,3,  and  then  assembles  them  into  a  subdecomposition  according 
to  (B33)  .  Besides  that,  the  three  elements  are  stored  temporarily, 
to  be  assembled  in  the  main  line  later  on  into  a  3x9  subdecomposi¬ 
tion  and  written  into  a  file,  each  with  three  indices  determined  by 
the  nodal  counts  o,n,m  in  a  formula  of  the  type  (B31) .  They  are 
used  later  to  assemble  the  matrix  of  the  advective  terms 

(WV  + WV  )  in  Equations  (5.18),  (5.19). 
onm  oimr  n 


>, 


3 . 9  Calculation  of  the  Total  Column  Number  in  the  Nonlinear  Matrix 

wv 

onm 

V 

The  column  number  N  ,  say,  of  the  nonlinear  matrix  W 

c  J  onm 

after  subdecomposition  can  be  calculated  in  the  following  way.  First, 
there  are  two  different  types  of  interaction  of  nodes  to  distinguish: 
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1.  self  interaction:  this  relates  to  matrix  elements  like  W  , 

onn 

where  the  node  r  interacts  with  itself.  As  has  been  shown  by 

n  J 

(B33)  ,  this  leads  to  a  six  column  subdecomposition  or,  equivalently, 
to  six  interactions  between  components. 

2.  mutual  interaction:  this  relates  to  elements  like  WV  ,  n  ^  m. 

onm 

If  r  and  r  are  on  different  tetrahedrons  only,  WV 

n  m  onm 

vanishes,  as  is  evident  by  the  disjoint  support  areas  of  the  basis 

functions  W  ,  W  .  This  means  no  interaction.  If  r  and  r 
r  r  n  m 

n  m 

are  on  one  and  the  same  tetrahedron,  this  evidently  leads  to  a 
nine  column  subdecomposition  or  to  nine  interactions  between 
components . 

In  a  discretized  rectangular  domain,  four  types  of  tetrahedrons 
are  possible. 

1.  Interior  tetrahedrons  with  four  counted  nodes. 

2.  Tetrahedrons  with  one  face  in  a  surface  plane  of  the  domain.  Three 
nodes  are  counted. 

3.  Tetrahedrons  with  two  faces  in  surface  planes,  i.e.,  tetrahedrons 
sitting  on  the  edges  of  the  domain.  Two  nodes  are  counted. 

4.  Tetrahedrons  with  three  faces  in  surface  planes,  i.e.,  tetrahedrons 
sitting  in  corners.  Only  one  node  is  counted. 

Clearly,  the  following  interactions  can  be  determined  for  these 


tetrahedrons : 
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Counted  nodes 

Intc 

Self 

tractions 

Mutual 

Total  Component 

1 

1 

0 

6 

2 

2 

1 

2x6+9  =21 

3 

3 

3 

3x6+3x9=45 

4 

4 

6 

4x6+6x9=78 

Moreover,  the  self  interaction  of  a  node  shared  by  two 
tetrahedrons  can  only  be  counted  once. 

This  clearly  gives  a  formula  for  the  total  number  of  inter¬ 
actions  in  a  given  discretized  domain  or,  equivalently,  the  total 

number  of  columns  N  ,  of  its  nonlinear  matrix  W^7  : 

c  onm 

Nc  =  T  •  6  +  T2  •  21  +  T3  •  45  +  T4  •  78  -  1*6  (B34) 

=  number  of  tetrahedrons  with  i  counted  nodes,  i  =  1,2, 3, 4 
1^  =  number  of  shared  faces 

3.10  Assembly  of  the  Nonlinear  Matrix  ^onm  in  the  Main  Line  Program 

The  main  line  consists  basically  of  five  nested  loops  (see 
Figure  15)  indexed  by  K,J,I,T,L1  in  the  order  from  outer  loop  to 
inner  loop.  K,J,I  iterate  through  the  discretized  values  of  the 
Z,Y,Z-axes  to  generate  the  rectangular  blocks,  which  split  into  five 
tetrahdrons  each.  T  goes  through  these  tetrahedrons  and  LI  through 
the  four  nodes  per  tetrahedron.  Clearly,  most  of  the  operations  have 
to  be  done  within  LI  (see  Figure  16) ,  as  is  the  case  for  the  assembly 


of  WV  .  This  is  done  according  to  the  structure  of 

onm  onm 


as  shown 
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in  Figure  14.  The  first  6N  columns  of  the  3N  x  Nc  matrix  are 
(partially)  occupied  by  elements  from  self-interacting  nodes  only, 
and  the  remaining  (N  -  6N)  columns  by  the  mutually  interacting  nodes. 
The  sparseness  of  this  matrix  (as  of  the  linear  matrices)  clearly 
depends  on  the  level  of  discretization;  the  finer  the  discretization, 
the  more  tetrahedrons  are  used  and  therefore  the  sparser  the  matrix 
will  be.  Since  the  logic  of  the  program  can  be  followed  by  the  diagrams 
of  Figures  15  and  16,  only  a  list  of  the  FORTRAN  symbols  used  in  the 
program  and  in  Figures  15,  16  as  well  is  given  here. 


Arrays : 

C(5,4),  D(5,4),  E(5,4)  -  X,Y,Z  -  coordinates  for  the  vertices  1  to  4 

of  the  tetrahedrons  1  to  5  per  block. 
SX(5),  SY(5),  SZ(5)  -  sums  of  the  X,Y,Z-coordinates  over  the  four 

vertices  of  the  tetrahedrons  1  to  5  per 
block . 

DET(5)  -  functional  determinants  between  barycentric 

and  cartesian  coordinates  for  the  tetrahedrons 
1  to  5  per  block. 

XNC(4,4),  YNC(4,4),  ZNC(4,4)  -  X,Y,Z  -  coordinates  of  the  centroids 

opposite  to  the  faces  1  to  4  per 
tetrahedron  1  to  4  per  block. 

TNL(3)  -  array  for  temporary  storage  of  the  three 

wr 

i  i  ^TTl 

elements  <W  W  — - —  >  ,  i  =  1,2,3  per 

r  i  r  1  X  .  r 

on  l 

nonlinear  subdecomposition. 

ANL(3,6,4,4)  -  array  for  temporary  storage  of  the  3x6  sub¬ 


decompositions  stemming  from  the  self  interactions . 
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ANLF (3N,N  ) 
ICl(Nc),IC2(Nc) 


NC(L3,T,I,J,K) 


NC0L(4,4) 


-  array  holding  the  3N  x  N  nonlinear  matrix  W  . 

J  &  c  onm 

-  arrays  assigning  the  indices  (IC1,IC2)  of  the 

appropriate  velocity  components  in  the  product 

vector  (V, . (r  ))  to  their  column  number  between 
li  n 

1  and  N  . 

c 

-  array  assigning  the  count  number  NC  of  a  node 
to  its  identity  (L3,T , I , J , K) . 

-  array  assigning  the  count  number  of  combinations 
of  a  pair  of  nodes  within  a  tetrahedron  to  the 
pair  of  their  identity  numbers  (1 .  .  .4, 1 .  .  .4) 
within  this  tetrahedrons.  Only  used  for 
combinations  on  mutual  interactions. 


Loops : 

I,J,K 

T 

L1,L2,L3 


LL2 , LL3 


-  increments  of  X, Y,Z-coordinates . 

-  tetrahedrons  1  to  5  per  block. 

-  3  nested  loops,  each  one  of  them  going  through 
the  four  nodes  per  tetrahedron. 

-  2  nested  loops  within  LI,  each  one  of  -them 
going  through  the  four  nodes  per  tetrahedron. 


3.11  Routines  ADV  and  TRNSP 

The  discretized  versions  (5.18),  (5.19)  for  second  and  third 
perturbation  differ  from  the  one  for  the  first  perturbation  by  the 
term: 


l 


r 

n 


,r  eS 
m 


(WV  +  WV  )V.  (r  )V.  (r  ) , 
onm  omn  l  m  j  n 


r  e  S 
o 


(B35) 
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with  indices  i,j  denoting  perturbations. 

Including  the  LHS-inverse  denoted  by  W  ^  and  introducing  fluid 

LH 

space  subdecomposition  as  well,  this  term  can  be  rewritten: 


n,m 


(WV  + WV  )V.  V. 
onm  omn  1m  jn 


(WV  ) 
on  m 


0  =  1 , . . . , 3N 


(B36) 


with  the  indices  o,n,m  going  from  1  to  3N  over  all  3N  velocity 
components  through  the  N  points  of  S. 

Moreover: 


(WV  ) 
on  m 


(WV  ) 

on  m 


onm 


omn 


■  3N  x  3N  matrix  indexed  by  m 


V  V 

(From  the  definition  (5.17)  of  W  it  is  clear  that:  (W  )  ^ 

v  onm  on;m 

V 

(W  )  .  Clearly,  the  two  sets  are  not  identical  since 

on  m  _  J 

WV  ^  wv  .) 

onm  omn 

This  reveals  the  term  in  the  square  brackets  of  (B36)  as  a  linear 

combination  of  matrices  and  W  ^  can  therefore  be  moved  into  the  sum: 

LH 


3N 

I 


n 


(WV  ) 
on  m 


V. 

jn 


(B37) 


V 

In  the  subroutine  ADV,  the  elements  of  the  subdecomposed  W  ,  as 

’  r  onm 

calculated  by  QMTRXN  and  stored  with  their  three  indices,  are  sorted 

once  after  the  index  n  and  once  after  the  index  m.  This  gives  all 

3N  subdecompositions  of  (W^  )  and  (WV  )  .  After  this  sorting 

r  on  m  on  m 
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V  V  - 1 

process,  the  pairs  (W  )  ,  (W  )  are  added  and  multiplied  by  WTI,  . 
r  r  on  m  on  m  r  LH 

A  transposing  routine  rearranges  the  elements  of  the  resulting  3N 
matrices  in  the  sequence  in  which  they  are  being  used  to  assemble  and 
evaluate  the  RHS  of  (B36)  in  the  programs  which  evaluate  Equations 

(5.18)  and  (5.19).  The  same  type  of  transposing  routine  (called 
TRNSP)  is  used  to  rearrange  the  elements  of  all  linear  matrices  into  the 
sequence  as  they  flow  into  routines  used  later  to  evaluate  Equation 

(5.19)  or  (5.20)  respectively. 


6N  Columns 
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3 
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cj 

o 
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r—  K  -  Z 


K 


J  -  Y  >  axis  is  incremented  by  discrete  values 

-  I  -  X 


Vertex  coordinates  C(5,4)D(5,4)E(5,4) ,  their 
sums  SX(5) ,SY(5) ,SZ (5) ,  functional  determinant 
DET(5)  and  nodal  coordinates  XNC(4,4),  YNC(4,4), 
ZNC(4,4)  are  determined  per  block  for  con¬ 
figuration  A  or  B. 


—  T  -  loop  iterates  through  the  5 
tetrahedrons 

Integrals  /T  XE'XYEyZEzdXdYdZ  =  QMLTI1  • 

DET(T)  are  calculated  for  Ex  +  Ey  +  Ez  <_  6 
and  stored  in  QT(7, 7,71.  SORT  2  and 
MATINV  calculate  coefficients  for  basis 
factors . 


-  LI  -  first  nodal  loop 

See  Figure  16  for  details. 


LHS-matrix  is  inverted  over  into  RHS-matrices  by  IMSLLIB  routines 
LINV3F  and  VMULFF . 


Figure  15.  Outer  Loops  of  DSCRB 
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LI 


skip  nodes  in  3D 


L2 


skip  nodes  in  3D 


L3 


skip  nodes  in  3D.  If  so,  go 
to  call  for  QMTRXA ,  QMTRXB 


count  nodes  via  the  first  encountered  set  of 
iterations  in  L3.  Check  node  identities  for 
previous  occurrences  and  skip  if  necessary. 


call  QMTRXN (LI , L2 , L3)  and  assemble  ANL(3,6),  TNL(3) 


count  mutual  interactions  of  nodes  via  first 
encountered  set  of  iterations  in  L2,L3  and  store 
in  NCOL (L2 , L3) . 


if  L3  =  4  call  QMTRXA ( . . . LI , L2 . . . ) , QMTRXB ( . . . LI , L2) 


L3 

L2 


copy  ANL(3,6)  into  part  I  of  ANLF(3N,Nc). 
Calculate  IC1,IC2  for  part  I  of  ANLF. 


LL2 


skip  nodes  in  3D 


LL3 


skip  nodes  in  3D 


redetermine  mutual  interactions  and  pick  their  counts 
from  NCOL.  Use  them  to  calculate  column  no.  in 
part  II  of  ANLF  and  distribute  TNL  accordingly. 
Calculate  IC1,IC2  for  part  II  of  ANLF. 


LL3 

LL2 

LI 


Figure  16.  Inner  Loops  of  DSCRB 
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1.  FLOWS  TESTED  FOR  BIFURCATIONS 

All  flows  show  the  main  characteristics  already  mentioned  in 
the  text:  linear  shear  flows  or  one  component  flows  undergo  stationary 
bifurcations  only;  nc  ^  |  v  |  ,  w0  ^  |  v  |  ,  ^  I  v  ^  I  »  ^2  ^  I  v  ^  I  > 

qualitatively  correct  behaviour  is  found  in  almost  all  cases  by  loss  of 
stability  with  increasing  shear.  All  flows  were  unstable  for  positive 
viscosities  except  flow  Number  3)  which  became  critical  for  v  =  1  at: 


U  =  5  ,  U  =  98 . 46y  , 
X  z 

for  v  =  .1  at: 

U  =  5.,  =  48 . 98y 

X  z 

and  for  v  =  .01  at: 


U  =  5.,  U  =  44 . 99y 

x  z 

U  =1.,  U  =  9 . 35y  . 

X  z 


In  all  four  cases  the  shear  coefficient  of  U  was  used  as  bifurcation 

z 

parameter  and  the  stationary  bifurcations  went  from  unstable  to  stable 
with  increasing  shear. 


Flows  with  Linear  Shear 


1)  U 

y 

2)  U 

x 

3)  U 

x 

4)  U 

x 


C  X 
yx 


C  ,  U  =  C  X 
x  y  yx 


C  ,  U  =  C  Y 
x  z  zy 


C  X  ,  U  =  C  Y  , 
xx  y  yy 
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5)  U 

x 


C  X 
xx 


U 

y 


c  Y 
yy 


u  =  c  z 

z  zz 


Only  stationary  bifurcations  could  be  found,  which  caused  a 

loss  of  stability  by  increasing  the  shear  coefficients  above  their 

critical  value.  These  occur  e.g.  for  viscosity  v  =  -1  and  for 

Flow  1)  at  C  =  6.56,  Flow  2)  at  C  =  .3,  C  =  7.1  and  Flow  3) 
yx  x  yx 

at  C  =  10,  C  =  52.8  .  The  critical  values  of  the  shear  coefficients 
x  zy 

show  only  a  weak  dependence  on  the  shear- free  components  such  as  C 

X 

in  Flow  2)  and  3)  and  in  Flow  4) . 

For  the  Flows  4)  and  5)  the  shear  coefficients  were  set  as 
linear  functions  of  a  bifurcation  parameter  n  and  the  occurrence 
and  value  of  critical  shear  coefficients  varied  with  the  combination  of 
proportionality  constants  used.  The  critical  shear  coefficients  however 
stayed  within  the  orders  of  ma°nii'ud^  of  10  or  1. 


One-Component  Flows  with  Nonlinear  Shear 

6)  U  =  C  Y2 

x  xy 

7)  U  =  C  Z3 

x  xz 

Again  only  stationary  bifurcations  could  be  found.  They  occurred 

for  viscosity  v  =  -1  and  Flow  6)  at  C  =  8.07  and  Flow  7)  at 
7  xy 

C  =9.39. 
xz 


Two-Component  Flows  with  Nonlinear  Shear 


8) 

9) 

10) 


U  =  C  Y  ,  U 
x  xy  z 


U  =  C  Z  ,  U 
x  xz  y 


U  =  C  -  C  Y  , 
x  x  xy 


C 

z 
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All  three  flows  show  unstable  Hopf  bifucations  (besides  station¬ 
ary  ones')  for  shear  coefficients  between  zero  and  one  for  flows  9)  and 
10)  and  around  five  for  flow  8)  (Viscosity  v  =  -1.)  They  are 
subcritical  for  flow  9) ,  supercritical  for  flow  10)  and  both  for 
flow  8) .  Flow  8)  also  showed  the  first  stable  (subcritical)  Hopf 

bifurcation  found  around:  v  =  -2.,  C  =  -11.,  C  =9.. 

z  xy 


Three  Component  Flows  with  Nonlinear  Shear 


11) 

U 

— 

2 

C  Y  , 

U  =  C  X, 

U  =  C 

X 

xy 

y  yx 

z  z 

12) 

U 

— 

2 

C  Y  , 

u  =  c  x3. 

U  =  C 

X 

xy 

y  yx 

z  z 

13) 

U 

— 

C  Z, 

U  =  c  x3, 

u  =  c  r* 

X 

xz 

y  yx 

z  zy 

14) 

U 

_ 

c  z3. 

U  =  c  x3. 

U  =  C  Y3 

X 

xz 

y  y * 

z  zy 

15) 

U 

— 

2 

C  Y  , 

u  =  c  z2. 

u  =  c  x2 

X 

xy 

y  xz 

z  zx 

16) 

U 

_ 

4 

C  Z  , 

4 

U  =  C  X  , 

4  * 

U  =  C  Y 

X 

xz  * 

y  yx 

z  zy 

17) 

U 

_ 

c  z 3  + 

C  Y3,  U  = 

3  3 

C  X  +  C  Z  , 

U 

=  C  YJ 

3 

+  C  X 

X 

xz 

xy  y 

yx  yz 

z 

zy 

zx 

18) 

U 

— 

2 

C  Y  + 

c  z3,  U  = 

2  2 

C  X  +  C  Z  , 

U 

2 

=  C  X 

3 

+  C  Y 

X 

xy 

xz  y 

yx  yz 

z 

zx 

zy 

For  all  flows  Hopf  bifurcations  were  quite  frequent.  In  general, 
the  critical  values  for  the  coefficients  stayed  in  their  order  of 
magnitude  between  one  and  ten,  irrespective  whether  some  of  them  were 
kept  fixed  (bifurcations  on  flows  with  some  of  their  coefficients  held 
fixed  at  other  orders  of  magnitude  were  rather  soarse)  or  all  were  linear 
functions  of  a  bifurcation  parameter  n  .  Also  secondary  bifurcations 

•k 

An  experimental  version  of  the  discretization  program  extended  to 
spatial  dependencies  up  to  power  6  was  used  here. 
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for  a  time  setting  of  t  =  0  on  the  first  orbit  were  found.  All 
possible  types  and  combinations  of  criticality  and  stability  could  be 
observed.  However,  supercritical  and  stable  bifurcations  appeared  more 
frequent  than  others. 


Flows  with  Affine  and  Power  Dependencies  on  the  Bifurcation  Parameter  n 


19)  Ux 
2°)  ux 
u 

z 

21)  U 

x 

22)  U 

x 

U 

y 

u 

z 

23)  U 

x 

U 

y 

u 

z 


(C  -n+CJY2,  u  =  (C  •  n  +  C„)z2 *,  u  = 
xy  5  y  yz  2  z 

el  2  e2  2 

(C  •  n  +  CJY  ,  U  =  (C  •  n  +  CL)Z  , 
xy  1  1J  y  K  yz  '  2  ’ 


(Czx'n  +  C  3)X 


(Czx-n  3  +  C3)x2 


n  (C  +  C  Y2) ,  U  =  n(C  +  C  •  Z2) , 
x  xy  y  y  yz 

el  2  e2  3 

(C  n  +CJYZ  +  (C  n  + 

xy  1  xz  2 

S3  3  e4  2 

(C  •  n  +  C-)X°  +  (C  n  +  CJZ 

^  yx  l  ^  yz  4 

-7 

(C  •  n  +  Cc)X  +  (C  n  +  C,)Y'5 

^  zx  1  5  zy  6 


U  =  n(C  +  C  •  X  ) 
z  z  zx  J 


C,  +  C  n  +  (C  +  C  n)Y' 

1  x  4  xy 

C„  +  C  n  +  (Cc  +  C  n)Z“ 

2  y  5  yz 


C  +  C  n  +  (C,  +  C  n)X 

3  z *  1  6  zx  J 


The  findings  for  this  group  are  in  principle  the  same  as  for 
the  previous  one  with  linear  dependencies  only.  However,  stable  and 
supercritical  Hopf  bifurcations  were  more  frequent  and  less  sensitive 
against  variation  of  parameters,  coefficients,  and  exponents  in  the 
affine  and  power  dependencies.  Besides,  third  level  bifurcations 

(time  settings  on  first  and  secondary  orbit  were  t  =  0)  could  be 
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found  only  in  this  group. 


APPENDIX  D 


1.  INTRODUCTION  TO  NUMERICAL  WORK  ON  THE  NS -EQUATION 

Most  of  the  numerical  work  done  so  far  on  the  NS-equations  is 
restricted  to  two  dimensional  flows,  since  there  the  NS-equations  can 
be  replaced  by  a  so-called  stream  function-vorticity  formulation. 

Applying  the  curl  operator  on  the  NS-equation  gives  the 
vorticity  equation: 


+  U*Voo  +  go‘VU  =  vV  oj, 

O  L 


(A) 


(with:  03  =  V  x  u  as  vorticity) 


(Al) 


and  eliminates  the  pressure  gradient  since  V  x  v  e  0.  In  a  two 
dimensional  flow  field  the  vorticity  vector  reduces  to  one  component, 
to  ,  say  perpendicular  to  the  plane  spanned  by  X,Y  of  the  flow  and 

Li 

the  vorticity  equation  can  be  written  as  a  scalar  equation: 


8u)„ 


at 


3oj, 


ao). 


=  -u 


x  ax 


u 


y  ar 


+  v 


2  2 

a  oj„  a  037 

Zj  Li 

— =—  +  — — 


(B) 


ax 


3Y 


au. 


with : 


03,  = 


Y 


ax 


5c 

aY 


(Bl) 


Also,  in  two  dimensions,  the  incompressibility  condition: 


div  U  =  0 


(C) 


reduces  to: 


au 

x 

ax 


au 


3Y 


0, 


(Cl) 
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which  can  always  be  satisfied  by  introducing  the  (scalar)  stream 
function  ip  as: 


136 


U 

x 


dip 

W  ’ 


u 

y 


dip 

dX 


(D) 


Clearly,  in  a  stationary  flow  the  following  relation  holds  along  the 
particle  trajectories: 


=  %  =>  -U  dX  +  U  dY  =  0 
U  dY  y  x 

y  3 


(E) 


and  the  isolines  of  the  stream  function  (i.e.  the  streamlines)  must  be 
identical  to  the  trajectories  since 


d*  -  If  dx+  §y  dY=  -Vx+  VY 


(F) 


Combining  the  two-dimensional  vorticity  (Bl)  with  the  stream  function 
(D)  gives  the  elliptic  Poisson  equation: 


2  2 

LA  +  LA  = 

3X2  9 Y2 


-co. 


(G) 


Equations  (B) ,  (D) ,  (G)  are  four  equations  for  the  four  unknowns 

U  ,  U  ,  \p ,  a)  .  They  are  mostly  solved  in  engineering  applications 
x  y  z 

under  the  appropriate  boundary  conditions  by  finite-difference  or 
finite  element  methods.  Their  advantages  as  well  as  their  disadvantages 
are  evident;  pressure  is  eliminated  and  the  incompressibility 
condition  is  incorporated  most  conveniently  steady-state  flow 
patterns  can  be  obtained  according  to  relations  (E) ,  (F)  most  easily 
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by  graphing  the  isolines  of  the  stream  function  ip  . 

On  the  other  hand,  however,  the  stream  function  concept  cannot 

be  applied  in  three  dimensional  flows  and  time- dependent  solutions  would 

require  substantial  modifications  to  introduce  evolution  equations  for 

U  ,  U  ,  and  \p  as  well.  Since,  in  the  present  context,  a  monitoring 
x  y 

of  the  evolution  time  of  a  velocity  vector  through  a  particular  point 
in  the  fluid  domain  is  required,  this  vorticity-stream  function 
concept  cannot  be  applied. 

The  publications  on  the  stream  function-vorticity  approach  are 
numerous.  See  e.g.  Morris  [15],  Smith,  Jr.  and  Kidd  [16],  Rubin  and 
Graves,  Jr.  [17,18],  Boney  [19],  Hefner  [20],  Hirsh  [21],  Suttles  [22], 
Zoby  [23]  for  treatment  and  description  of  various  technical  and 
computational  details,  the  associated  problems  of  convergence  and 
the  numerical  algorithms  applied  to  achieve  convergence.  Most  of  the 
work  deals  with  the  so-called  cavity  problem,  where  the  fluid  is 
considered  enclosed  in  a  rectangular  cavity  and  its  steady  motion  is 
driven  by  one  sliding  wall  or  more  via  non  slip  boundary  conditions. 

In  the  seventies,  finite-element  techniques  in  connection  with 
the  Galerkin  method  were  adapted  to  solve  the  NS-equations .  However, 
the  amount  of  numerical  operations  and  storage  requirements  are 
considerably  higher  than  in  e.g.  finite-difference  methods,  since  the 
coefficient  matrices  contain  more  nonzero  entries.  Representative 
papers  would  be,  e.g.,  Cheng  [24],  Taylor  and  Hood  [25],  Kawahara  et.  al. 
[26],  Fortin  and  Thomasset  [27],  Olson  and  Tuann  [28]. 

Again  most  of  the  work  done  so  far  is  restricted  to  two 
dimensions.  One  of  the  reasons  for  this  seems  to  be  that  the  amount  of 
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numerical  computations  and  storage  requirements  increases  by  almost 
one  order  of  magnitude  from  two  to  three  dimensions,  as  could  be 
concluded  in  the  present  work.  Moreover,  the  finite-element  technique 
depends  on  a  triangulation  (i.e.,  partition  into  triangles)  of  the  two- 
dimensional  fluid  domain,  which  naturally  has  to  be  based  on  tetra¬ 
hedrons  in  three  dimensions  instead  of  triangles.  The  problems 
involved  therewith  do  not  permit  an  almost  straight-forward  inter¬ 
polation  from  two  to  three  dimensions,  as  is  the  case  for  rectangular 
finite  differences.  However,  application  of  the  finite-element/Galerkin 
technique  to  three  dimensions  seems  at  present  the  only  feasible 
method  to  discretize  the  NS-equation  for  the  problem  at  hand  into  a 
system  of  evolution  ODE's  which  can  be  solved  without  algorithms  for 
correction  to  satisfy  incompressibility  in  between  two  timesteps,  as 
is  important  for  the  applicability  of  the  program  BIFOR  2  ,  which  tests 
the  ODE-system  for  Hopf  bifurcations.  The  main  reference  for  the 
analytic  background  on  which  the  here  employed  techniques  are  based,  as 
well  as  for  a  very  thorough  and  complete  extensive  analysis  of  the  NS- 
equations  is  the  book  by  Temam  [12] . 


APPENDIX  E 


1.  THE  INERTIAL  SUBRANGE  LAW 

(a)  Kolmogorov’s  Derivation  [13],  1941 

Kolmogorov  (13)  divided  the  spectrum  of  energy  per  unit  mass 
versus  wave  number  of  an  isotropic  homogeneous  turbulent  flow  into 
three  ranges: 

1.  the  production  range  at  low  wave  numbers,  which  draws  energy  from 
the  mean  flow  and  converts  it  into  turbulent  fluctuations; 

2.  the  inertial  range  at  intermediate  wave  numbers,  in  which  energy 
is  shifted  upwards  towards  higher  wave  numbers  without  much 
production  or  dissipation  taking  place; 

3.  the  dissipation  range  at  high  wave  numbers,  in  which  the  energy  put 
into  the  production  range  is  finally  dissipated  into  a  motion 
comparable  to  Brownian  motion. 

Assuming  that  the  transfer  of  energy  from  production  range 

through  inertial  range  into  dissipation  range  is  a  local  phenomenon, 

these  three  ranges  are  virtually  independent  of  each  other.  In 

particular,  the  spectra  in  the  inertial  and  dissipation  ranges  will  not 

be  much  affected  by  the  precise  way  in  which  energy  is  put  into  the 

production  range.  Moreover,  in  a  stationary  mean  flow  the  production 

rate  of  turbulent  energy  must  be  equal  to  the  dissipation  rate,  e  ,  say. 

Calling  the  spectral  density  E(k)  and  recognizing  it  as 

energy  per  unit  mass  within  the  wave  number  range  (k,k  +  dk),  it  must 

3-2 

have  dimension  L'  T  .  The  dissipation  rate  e ,  as  energy  per  unit  mass 

2  -3 

per  time,  must  have  dimension  L  T  .  It  is  also  clear  that  smaller 
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eddies  (i.e.  larger  k)  contain  less  energy  in  their  associated 
wave  number  range.  One  can  therefore  write: 

E(k)  =  Ck  •  fk(k)  +  k,  '  (El) 


i.e.,  E(k)  decreases  monotonically  with  increasing 
On  the  other  hand,  a  higher  dissipation  rate 
more  energy  has  to  pass  through  from  k  to  k  +  dk. 


wave  number  k. 
e  means  that 
That  means : 


E  (k)  =  C  •  f  (€)  +  e. 


(E2) 


i.e.,  E(k)  increases  monotonically  with  increasing  dissipation  rate 

€  . 

It  is  also  evident  that  in  the  limit  for  laminar  flow  €  ->  0  : 

E(k)  =  0  for  e=0  (E3) 


Since  dissipation  in  the  inertial  subrange  is  negligible,  E(k)  cannot 
depend  on  the  viscosity  there.  Therefore  the  only  possible  dependent 
parameters  left  are  6  and  k.  E(k)  must  therefore  be  a  combination 
of  (El)  and  (E2)  under  consideration  of  (E3) ,  which  yields :. 

E (k)  =  C  •  f£(e)  •  fk(k)  (E4) 


the  functions  f  (e)  and  f  (k)  can  now  be  easily  dimensionally 

6  K. 

scaled  as: 


f 

6 


(e) 


fk(k)  =  k 


-5/3 


(E5) 


and  C  as  a  dimensionless  constant  to  obtain  the  correct  dimension 


V. 
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3-2  2-3  -1 

L  T  for  E(k)  and  observing  L  T  and  L  as  dimensions  for  e 

and  k  respectively. 

(E4)  and  (E5)  give  now  the  well-known  inertial  subrange  or 
(-5/3) -law: 


E(k)  =  C  •  e2/V5/3  (E6) 

b)  Heisenberg’s  Derivation  [14],  1948 

This  derivation  is  directly  related  to  the  finite  closure  methods. 
Beginning  with  the  NS-equation 


_9_ 
9 1 


U.(X) 


[U.cxpu.cxp]  =  vV^U.QCp 


apCxp 


(E7)  * 


where  the  incompressibility  condition 


7—  u  (X  )  =  0  (E8) 

lj  3 


is  used  in  the  advective  term,  the  equation  (E7)  is  multiplied  with 
^k^2^  an<^  ensemble  averages  are  taken: 


fi  WW  +  iiq-  WWW 


(E9) 


•  vVi  WW  -  sqr  WW 


Since  homogeneity  is  assumed,  the  averaged  products  depend  on  the  spatial 
difference  r  e  X9  -  X^  only.  Defining  the  correlations: 


For  reasons  of  conciseness,  time  dependence  is  not  written  explicitely. 
Subscripted  numbers  denote  points  and  letters  denote  components. 
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R. . (r)  e  U. (X,)U.X_ 
ij  1  1'  j  2 


S.  (0,r)  =  S.  (r)  =  U.  (XJU.  (XJU.  (X0) 

ijkv  ljk  l  1  j  1  kv  2J 


p.(r)  =  P(x1)u.  (x2; 


(E10) 

(Ell) 

(E12) 


and  using 


3r. 


3X 


,  Equation  (E9)  writes  as 


li 


h  Rik(?)  -  w:  sijk(7)  =  vVikc7)  *  a!:  V7) 

J  r  i 


(E13) 


In  the  same  way,  the  following  conditions  can  be  derived  straight 
forwardly  from  the  incompressibility  condition  (E8) : 


3  R..  (?)  =  4k  s.  (?)  -  4k  M?)  =  0 


3r.  ik 

l 


3rk  ijk 


3r .  i 
i 


(E14) 


Fourier-transforming  the  correlations  (E10) ,  (Ell),  (E12) 


qyio  =  (2T)'3  /  R.  ,(?)e"lkrd  Vol 
Y  (k)  =  C2-7T)"3  /  S  (?)e'lkrd  Vol 


0.00  =  (2tt)'3  /  P.(?)e'lkrd  Vol, 


(E15) 

(E16) 

(E17) 


gives  for  Equations  (E13)  ,  (E14) : 


£*.k(k)-ikjY..kCk)  =  -vk2*.k(k)  ♦  ik.ek(k) 


(E18) 


k  .  d) . ,  =  k,  y.  =  k.e.  =  0 

iTik  k'ljk  i  i 


(E19) 


Besides  homogeneity,  another  important  simplification  is  isotropy.  The 
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correlation  tensors  d> .  . ,  y.  ,  0. 

ij  ljk  1 


therefore  take  the  following  form: 


<t>.k(k)  =  a1(k)«.k  *  a2(k)kkkk 


y..v  (k)  =  b1(k)k.k.k  +  b  (k)k.6.,  +  b_(k)k.<S..  +  b  (k)6.. 
ijk^  lv  i  j  k  2  ;  l  jk  3W  j  lk  J  ij 


(E20) 
(E21 ) 


In  the  tensor  context  isotropy  is  synonymous  with  invariance 
under  the  full  rotation  group,  ie.,: 


^ik  ^ik  ajiaLk^jL 


Y. 


=  Y  —  (JL  CM  C*.  V 

ijk  ijk  Li  mj  nk  Lmn  ' 


=  aT,am,a  ,Yl 


(E22) 

(E23) 


with  a^  as  the  matrix  performing  the  rotation  of  coordinates  and  the 
prime  denoting  a  rotated  coordinate  system.  The  isotropy  conditions 
(E22)  ,  (E23)  can  now  easily  be  verified  for  the  structures  (E20)  and 
(E21)  respectively,  by  using  linearity  of  the  rotation  and  isotropy 
of  6^  .  Also  the  following  symmetry  is  evident: 


S 


(E24) 


The  pressure  correlation 


e . 
j 


vanishes 


0j(k)  =  o. 


(E  2  5 ) 


as  is  evident  by  physical  arguments,  since  under  isotropy  pressure  cannot 
correlate  with  a  particular  direction  such  as  the  velocity  vector  U. 

The  incompressibility  conditions  (E19)  simplify  the  correlation 
tensors  <Jk  ^  ,  ^ijk  father,  since  they  require: 
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a  (k)  +  k  a2(k)  =  0 


k  bx(k)  +  2b2(k)  =  b4(k)  =  0, 


(E26) 

(E27) 


And  can  now  be  expressed  by  only  one  scalar  function  each: 


4) .  •  (k)  =  (k2(5-r  -  k .  k  ) 

it  .  .  2  ik  i  k 

J  4trk 


(E28) 


Yijk(k)  -  iY(10  k.k.kk-ik2(V.k  +  k.Sik) 


(E29) 


Entering  (E25) ,  (E28) ,  (E29)  into  (El 8)  gives 


E(k)  +  vk2E(k)  =  27rk6y(k) 

o  c 


(E30) 


Forming  the  trace  of  Equation  (E28)  leads  to: 


E(k)  =  27rk  <j>ii (k) 


(E31) 


On  the  other  hand,  the  kinetic  energy  density  of  the  flow  6j^n  is: 


e,  .  =  x  R..  (0)  =  ■§■  U.  (XJU.  (X.) 
km  2  ii v  2  iv  1J 


=  §  /  |  (k)eilrd\L  =  |  /  ♦ijWd^k  , 

r=0 


(E32) 


and  substituting  (E31)  into  (E32)  gives: 


£kin  '  /  E«dk 


(E33) 


E(k)  must  be  therefore  the  spectral  density  of  energy. 
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Next,  one  observes 


/  k . y . . . (k) d\  =  0 
J  3  iji 


(E34) 


This  can  easily  be  derived  from  the  Fourier- transform: 


■  r  i  ikr  ,3^  9  c 

J  j  iji  9r .  ij l v 

J 


(E35) 


3  U .  (X,  ) U .  (1 ) U.  (X  +t)  =  U .  (X,  ) U .  (X , }  U.  (X,  +?)  , 


9r.  i  1  j  1  i  1  i  1  j  1  9X.  .  i  1 

J  lj 


For  r  =  0  and  the  incompressibility  condition 


9U.  (jy 

9  X  1  . 
lj 


this  last  term  gives 


1  _JL 

2  9X, 


(U.cxpu.  (xpu.cxp) 


1  ^ 

2  9X,  . 

lj 


U.  (X_)U.  (XJU.  (XJ  =  0 

1  1  J  1  1  1 


(E36) 


since  the  differentiation  takes  place  on  a  constant  value. 

Integrating  Equation  (E30)  over  k-space  therefore  gives: 


CO  OO  CO 

j  E (k)dk  +  v  /  k2E(k)dk  =  2tt  /  k6Y(k)dk  =  0  (E37) 

dt  0  0  0 


Heisenberg’s  assumption  is  now  that  the  loss  of  kinetic  energy  of  small 
eddies  (i.e.,  large  wave  numbers)  is  negligible  compared  to  the  loss  on 
large  eddies  (i.e.,  small  wave  numbers).  That  is,  under  use  of  (E33) ; 


. 
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It  ekin  ^k=  large  up  to  °o)  =  pv  /  E(k?)dk'  =  0 

k 


(E38) 


Integrating  (E30)  from  0  to  k  and  from  k  to  °°  : 


k  k  k 

~  !  E (k ' ) dk  '  +  v  /  k'2E(k')dk’  =  2t t  /  k'6y(k')dk'  (E39) 


0 


0 


0 


jr  I  E  (k '  )dk '  +  v  /  k'2E(k')dk'  =  2tt  /  k'6y(k')dk' 

k  k  k 


(E40) 


Following  the  finite- closure  method,  the  third-order  correlation  term 
with  y(k)  in  (E39)  is  now  to  be  replaced  by  a  second-order 
correlation  term  with  E(k).  A  relation  between  these  two  terms  is 
given  by  (E40)  and  (E38) ,  however,  for  large  wave  numbers  only. 

Based  on  the  viscosity  term  of  (E39)  ,  a  turbulence  viscosity 
n(k)  can  be  introduced: 


uu  _  rv 

vf  k '  2E  (k ' )  dk  '  =  ri"(k)  /  k'2E(k')dk' 

k  0 


(E41) 


recalling  from  (E37)  : 


lx  ^ 

—  2tt  /  k^yCk^dk'  =  2tt  /  k,6y(k')dk', 

0  k 


(E42) 


this  relation  (E42)  together  with  (E41)  and  (E38)  can  be  used  in 
(E40)  to  obtain: 


1c  k. 

2tt  /  k,6y(k')dk'  =  n(k)  /  k  '  2E  (k  ' )  dk ' 

0  0 


(E43) 


V* 


a  relation  between  second  and  third-order  correlations  expressed  by 
E(k)  and  y(k),  respectively,  for  small  wave  numbers.  Equation  (E39) 
now  becomes : 
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/  E(k')dk'  =  -Cn(k)  +  v)  /  k,2E(k')dk'  (E44) 

3t  0  0 

From  the  kinetic  theory  of  gases,  the  following  proportionality 
between  viscosity  v  and  the  kinetic  energy  per  unit  mass  is  known: 

i, 

v  'v  (kmT)'2  (E45) 

with  kmT  as  kinetic  energy  per  unit  mass  determined  by:  k  =  Boltzmann 
constant,  m  =  molecular  mass,  T  =  temperature. 

This  relation  between  energy  and  viscosity  can  now  be  used  to 
express  the  turbulent  viscosity  n(k)  in  terms  of  spectral  energy 
density  E(k).  Establishing,  as  in  Kolmogorov's  derivation,  the 
dimensions  of  the  involved  quantities: 


spectral  energy 

density: 

E(k)  dimension: 

3  -2 
St 

-1 

wave  number: 

k 

s 

_ 

c2  -1 

viscosity: 

n 

S  t 

energy  per  unit 

mass : 

/  E  (k)dk 

2  -2 
St 

the  turbulent  viscosity  can  be  dimensionally  scaled  by  use  of  the 
proportionality  (E45)  as: 

oo 

n(k)  =  /  n(k')dk’ 
k 


(E46) 


- 
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n  00 


where  a  is  a  dimension  less  constant. 
Equation  (E44)  becomes  now: 


(E47) 


jr  I  E  (k T )  dk  '  e  e 
0 


E(k') 


h 

dk '  +  v 


/  k ' 2E (k ' ) dk ' 
0 


(E48) 


Assuming  time  intervals,  where  the  loss  of  energy  of  large  eddies,  i.e., 
e  ,  can  be  considered  constant,  this  equation  is  solved  by: 


E  (k) 


C(£)2/3k"5/3 


8  4 

1  +  (k/k  ) 

3a  3 


14/3 


(E49) 


with :  k .  =  . 

J  \  3 

\  v 


C  \h 


Therefore: 


E(k) 

.  -5/3 
^  k 

for 

k  << 

k. 

1 

E(k) 

,  -7 

'v  k 

for 

k  >> 

k. 

J 

Besides  the  inertial  subrange  law  for  k  <<  k_.  Heisenbergs  derivation 
gives  a  power  law  for  the  dissipative  range  k  >>  k.  as  well.  This 
is  possible,  due  to  the  use  of  kinetic  gas  theory  in  the  derivation. 


2.  THE  TAYLOR  HYPOTHESIS 

Measurements,  which  are  to  be  used  to  determine  the  spectrum 
E(k)  introduced  in  Part  1  must  consist  of  a  series  of  measuring  devices 
distributed  one-dimensional ly  and  equi-spaced  with  a  certain  distance  d 
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parallel  to  the  mean  steady  and  homogeneous  flow  U  .  Moreover,  this 
set  of  measuring  devices  should  drift  along  with  the  mean  flow  (and  mean 
velocity)  U  to  make  sure  only  the  turbulent  fluctuations  are  measured. 
Evidently,  such  a  set  of  devices,  which  is  large  enough  in  number  and 
is  fixed  in  a  Lagrangean  coordinate  system  moving  with  the  mean  flow  is 
virtually  impossible  to  realize  in  practice;  and  measurements  have  to 
be  made  on  the  basis  of  an  Eulerian  system  with  the  measuring  devices 
fixed  w.r.t.  the  ground.  This  Eulerian  frame  can  now  be  linked  to 
the  desired  Lagrangean  frame  on  the  basis  of  a  hypothesis,  originally 
due  to  G.I.  Taylor,  that  the  sequence  of  turbulent  fluctuations  at  a 
fixed  point  is  statistically  the  same  as  if  the  spatial  pattern  of 
velocities  were  suddenly  frozen  and  swept  past  the  measuring  device  with 
the  speed  U  of  the  mean  flow.  The  set  of  several  measuring  devices 
fixed  within  a  Lagrangean  frame  as  described  above  can  now  be  replaced 
by  only  one  device  in  the  Eulerian  frame  which  measures  in  equi-spaced 
time  intervals  t  which  relate  to  the  equi-spaced  distances  d  of 
the  Lagrangean  frame  in  the  following  obvious  way: 

d  =  Ut  ,  .  (E50) 

which  is  the  mathematical  formulation  of  the  Taylor  hypothesis. 

Clearly,  the  Taylor  hypothesis  is  important  for  the  present  work, 
where  only  a  rather  limited  number  of  points  (i.e.,  four)  is  monitored 
over  a  time  period  in  an  Eulerian  frame  fixed  in  fluid  space.  Therefore, 
a  Fourier  analysis  can  only  be  done  w.r.t.  time,  which  gives  a 
frequency  spectrum.  However,  as  is  evident  from  (E50) ,  the  Taylor 
hypothesis  holds  in  the  context  of  wavespace  and  frequency  as  well: 
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k  =  f/U 


(E51) 


with  k  and  f  as  wave  number  and  frequency  respectively. 


APPENDIX  F 


1.  TRAJECTORY  PROJECTIONS  AND  SPECTRAL  GRAPHS 


O 


denotes  zero  on  trajectory  scales 
denotes  beginning  of  trajectory 
denotes  end  of  trajectory 


No  spectral  graphs  SI,  S2,  S10,  Sll  are  provided, 
since  no  spectral  analysis  was  performed  for  the  relating 
trajectories  of  Figures  T1 ,  T2,  T10,  Til.  All  trajectories 
belong  to  Flow  IV  except  T1  and  T2,  which  belong  to  Flow  III. 
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FIGURE  T1 

Bifurcation  parameter  n  .55 

No.  of  timesteps  851 

Stepsize  .03 

Init.  condition  for  .001  (all) 

Average  size  of  spatial  random  perturbations  n/a 

Average  no.  of  timesteps  between  random  perturbations  n/a 

Period  of  1.66 
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FIGURE  T2 


Bifurcation  parameter  r\  .55 

No.  of  timesteps  2000 

Stepsize  .01 

Init.  condition  for  V9  .001  (all) 

Average  size  of  spatial  random  perturbations  n/a 

Average  no.  of  timesteps  between  random  perturbations  n/a 

Period  of  1.66 
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FIGURE  T3 


Bifurcation  parameter  n  .3 

No.  of  timesteps  1000 

Stepsize  .03 

Init.  condition  for  .001  (all) 

Average  size  of  spatial  random  perturbations  n/a 

Average  no.  of  timestepts  between  random  perturbations  n/a 

Period  of  1.68 

Mean  Velocity: 
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FIGURE  T4 
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FIGURE  T5 


Bifurcation  parameter  n  .325 

No.  of  timesteps  500 

Stepsize  .03 

Init.  condition  for  V2  .001  (all) 

Average  size  of  spatial  random  perturbations  n/a 

Average  no.  of  timesteps  between  random  perturbations  n/a 

Period  of  V  1.66 

Mean  velocity: 

Point  1  1.95 
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FIGURE  T6 


Bifurcation  parameter  n 
No.  of  timesteps 
Stepsize 

Init.  condition  for 

Average  size  of  spatial  random  perturbations 
Average  no.  of  timesteps  between  random  perturbations 
Period  of 
Mean  velocity: 
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FIGURE  T7 


Bifurcation  parameter  n 
No.  of  timesteps 
Stepsize 

Init.  condition  for 

Average  size  of  spatial  random  perturbations 
Average  no.  of  timesteps  between  random  perturbations 

V 

Period  of 
Mean  velocity: 
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FIGURE  T8 


Bifurcation  parameter  n 
No.  of  timesteps 
Stepsize 

Init.  condition  for  V9 

Average  size  of  spatial  random  perturbations 
Average  no.  of  timesteps  between  random  perturbations 
Period  of 
Mean  velocity: 
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FIGURE  T9 

Bifurcation  parameter  n  .47. 

No.  of  timesteps  1000 

Stepsize  .04 

Init.  condition  ror  V_  ^ 

2  rest  =  0 

Average  size  of  spatial  random  perturbations  n/a 

Average  no.  of  timesteps  between  random  perturbations  n/a 

Period  of  1.53 

Mean  velocity: 

Point  1  1.71 

2  1.76 

3  2.78 
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FIGURE  T10 


Bifurcation  parameter  n 

.32 

No.  of  timesteps 

500 

Stepsize 

.03 

Init.  condition  for 

.001  (X-comp .  Pt .  1) 
rest  =  0 

Average  size  of  spatial  random  perturbations 

n/a 

Average  no.  of  timesteps  between  random  perturbations 

n/ a 

Period  of 
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Mean  velocity: 
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FIGURE  Til 


Bifurcation  parameter  r\  .32 

No.  of  timesteps  379 

Stepsize  .03 

Init.  condition  for  V9  0  (all) 

Average  size  of  spatial  random  perturbations  .005 

Average  no.  of  timesteps  between  random  perturbations  2 

Period  of  1.66 

Mean  velocity: 

Point  1  1.95 

2  1.64 

3  2.65 
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FIGURE  T12 


Bifurcation  parameter  n  .27 

No.  of  timesteps  1000 

Stepsize  .04 

Init.  condition  for  .1  (all) 

Average  size  of  spatial  random  perturbations  .01 

Average  no.  of  timesteps  between  random  perturbations  1.5 

Period  of  1.7 

Mean  velocity: 

Point  1  1.68 

2  1.56 

3  2.43 

4  1.56 


(First  233  timesteps  are  unperturbed.) 
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FIGURE  T13 


Bifurcation  parameter  ri  „  28 

No.  of  timesteps  1000 

Stepsize  .04 

Init.  condition  for  V?  .01  (all) 

Average  size  of  spatial  random  perturbations  .01 

Average  no.  of  timesteps  between  random  perturbations  2 

Period  of  1.69 

Mean  velocity: 

Point  1  1.68 

2  1.55 

3  2.42 

4  1.56 


(First  233  timesteps  are  unperturbed.) 
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FIGURE  T14 

Bifurcation  parameter  n  .3 

No.  of  timesteps  500 

Stepsize  .03 

Init.  condition  for  V ^  0(all) 

Average  size  of  spatial  random  perturbations  .02 

Average  no.  of  timesteps  between  random  perturbations  2 

Period  of  1.68 

Mean  velocity: 

Point  1  2.20 

2  2.03 

3  2.95 

4  1.85 
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FIGURE  T15 


Bifurcation  parameter  n  ,325 

No.  of  timesteps  500 

Stepsize  .03 

__ 

Init.  condition  for  V  0  (all) 

Average  size  of  spatial  random  perturbations  .015 

Average  no.  of  timesteps  between  random  perturbations  2.5 

Period  of  1 . 66 

Mean  velocity: 

Point  1  1.95 

2  1.64 

3  2.65 

4  1.60 


YCOMP  .  POINT  2  YCOMP  .  POJNT  1 

.  jO  -O.IS  O.Ol  0.17  -0.90  -0.49  -0.07  0.34 


T15 (1) 


251 


-C  . b2 

1 - 


-C  .09 

i — e- 


C  .  34 


c .  n 


i  •  C  1 

f 


o 


4>« 

o 


o 
+  ' 


4 - 1 - $ - 1 - I - h 

-0 .92  -0  .09  0-34  C  .77  1.^1 

XC0MP  .  D0INT  1 


-  0 . 1 0  0  .  1 1 

4 - 4> - 1 - 


C  .  3 1  C.bl  0.72 

4 - 1 - r 


<> 


o 

I 


o 

o 


•> 


o 

i 


CD 


O 


4 - 0 - 1 - 1 - 1 - 

- 0 . 1 0  0  ■  1 1  c  .  3 1  C.bl 

XCCMP  .  D0!NT  2 


f 

C  .  )2 


GT0-  Ob’ 


115(2) 


252 


^ - 1 - 

-G  .7  4  -C.47 

XCOMP.  °0  I  NT  3 


H - $ - 1 - H 

C .20  G ■ C6  C  .33 


-C  -81 


-C  .55 


-C  .29 


+- 


-t- 


-C  . C3 

-t— $ - 


CO 

o 


o 

u 


o 

I 


T15 (3) 


253 


4 - 1 - 0 - 1 - 1 - 

-0-52  -0  .09  C  •  3  4  C-77 

XGGMP  .  p 0 1 N T  1 


in 

C5 


-0  .  i0 
H - 


o .  n 


0.31 

H - 


C  .  b  1 


C  . 


JO 

o 


o 

CL 


Q_ 

n 

o 

o 
~vj  -o 


o 

i 


CD 
.  _  I 


^ - $ - 1 - 1 - 1 - K 

-0.i0  0  .  1  1  C  .  3  i  c.bi  C./2 

XC0MP  •  °C  I  N T  2 


T15  (4) 


254 


4 - 1 - 

-C.74  - C  .47 

xconp  .  °0 I  NT  3 


H - - 1 - 1- 

-0.20  C.C6  C.j3 


T15 (5) 

-CS.90  -C.59  -G.Z8  C.C3  C.^4 

cn  4 - 1 - 1 - | - 1-  a ^ 


255 


zrorip  .  POINT  4  ZCOMP  -  POINT  3 

TO  0  .  C  4  0.18  0.32  0.46  -0.24  -0.i2  0.00  C  -  1 3  0.25 


T15  (6) 


i - 1 - <*> - 1 - 1 - h 

- C -19  -C.C3  0.12  C  .  28  C.43 

TCOMP  .  °OINT  3 


-0.36  -0.19  -C  .02  C  .  1 5  0.33 

1 - 1 - 1 — £ - 1 - r  to 


-n 


o 


<> 


o 


CO 

4-  1 


4 - 1 - 

-0.36  -G  -  i  9 

TC0MP  .  POINT  4 


+ 


C  -00 


- 

-0 .02 


1 - 

C  ■  *5 


257 


FIGURE  T16 

Bifurcation  parameter  n  .33 

No.  of  timesteps  645 

Stepsize  .03 

Init.  condition  for  V0  ^ 

2  rest  =  0 

Average  size  of  spatial  random  perturbations  .01 

Average  no.  of  timesteps  between  random  perturbations  2 

Period  of  1.65 

Mean  velocity: 

Point  1  1.75 

2  1.58 

3  2.43 

4  1.59 


First  200  timesteps  are  unperturbed. 
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FIGURE  T17 


Bifurcation  parameter  n  .445 

No.  of  timesteps  500 

Stepsize  .03 

Init.  condition  for  V0  .1  (all) 

Average  size  of  spatial  random  perturbations  .01 

Average  no.  of  timesteps  between  random  perturbations  2 

Period  of  1.55 

Mean  velocity: 

Point  1  1.61 

2  1.68 

3  2.67 

4  1.76 
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FIGURE  S3 


Bifurcation  parameter  n  .3 

No.  of  timesteps  1000 

Stepsize  .03 

Init.  condition  for  V  .001  (all) 

Average  size  of  spatial  random  perturbations  n/a 

Average  no.  of  timesteps  between  random  perturbations  n/a 

Period  of  V  1.68 

Mean  velocity 

Point  1  2.20 

2  2.03 

3  2.95 

4  1.85 


Sample  size/subsample  size  =  no.  of  subsamples 


1024/128  =8 


E./u.mss*wRVE#  E./u.mss*wnvE# 
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FIGURE  S4 


Bifurcation  parameter  n  .449 

No.  of  timesteps  800 

Stepsize  .04 

Init.  condition  for  V,  .001  (X-conrp.Pt .  1) 

2  rest  -  0 

Average  size  of  spatial  random  perturbations  n/a 

Average  no.  of  timesteps  between  random  perturbations  n/a 

Period  of  1.55 

Mean  velocity: 

Point  1  1.38 

2  1.76 

3  2.84 

4  1.95 


Sample  size/subsample  size  =  no.  of  subsamples 


896/128  =  7 
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FIGURE  S5 


Bifurcation  parameter  n  .325 

No.  of  timesteps  500 

Stepsize  .03 

Init.  condition  for  V0  .001  (all) 

Average  size  of  spatial  random  perturbations  n/a 

Average  no.  of  timesteps  between  random  perturbations  n/a 

Period  of  1.66 

Mean  velocity: 

Point  1  1.95 

2  1.64 

3  2.65 

4  1.60 

Sample  size/subsample  size  =  no.  of  subsamples  512/64  =  8 
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FIGURE  S6 


Bifurcation  parameter  n 
No.  of  timesteps 
Stepsize 

Init.  condition  for  V 

Average  size  of  spatial  random  perturbations 
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Mean  velocity: 
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FIGURE  S7 


Bifurcation  parameter  n 

A 

No.  of  timesteps 

800 

Stepsize 

.04 

Init.  condition  for 

(X-comp . Pt . 1) 
rest  =  0 

Average  size  of  spatial  random  perturbations 

n/a 

Average  no.  of  timesteps  between  random  perturbations 

n/ a 

Period  of  V 

1.60 

Mean  velocity: 
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2.61 

4 
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FIGURE  S8 


Bifurcation  parameter  n 
No.  of  timesteps 
Stepsize 

Init.  condition  for  V? 

Average  size  of  spatial  random  perturbations 
Average  no.  of  timesteps  between  random  perturbations 
Period  of  V 
Mean  velocity: 

Point  1 
2 
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4 

Sample  size/subsample  size  =  no.  of  subsamples 
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FIGURE  S9 


Bifurcation  parameter  n 
No.  of  timesteps 
Stepsize 

Init.  condition  for  V 

Average  size  of  spatial  random  perturbations 
Average  no.  of  timesteps  between  random  perturbations 
Period  of 
Mean  velocity: 

Point  1 
2 

3 

4 

Sample  size/subsample  size  =  no.  of  subsamples 
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FIGURE  SI 2 


Bifurcation  parameter  n  .27 

No.  of  timesteps  1000 

Stepsize  .04 

Init.  condition  for  V?  .1  (all) 

Average  size  of  spatial  random  perturbations  .01 

Average  no.  of  timesteps  between  random  perturbations  1.5 

Period  of  1.7 

Mean  velocity: 

Point  1  1.68 

2  1.56 

3  2.43 

4  1.56 

Sample  size/ subsample  size  =  no.  of  subsamples  768/64  =  12 


(First  233  timesteps  are  unperturbed.) 
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FIGURE  S13 


Bifurcation  parameter  n  .28 

No.  of  timesteps  1000 

Stepsize  .04 

Init.  condition  for  V  .01  (all) 

Average  size  of  spatial  random  perturbations  .01 

Average  no.  of  timesteps  between  random  perturbations  2 

Period  of  1.69 

Mean  velocity: 

Point  1  1.68 

2  1.55 

3  2.42 

4  1.56 

Sample  size/subsample  size  =  no.  of  subsamples  768/64  =  12 


(First  233  timesteps  are  unperturbed.) 
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FIGURE  S14 


Bifurcation  parameter  n  .3 

No.  of  timesteps  500 

Stepsize  .03 

Init.  condition  for  0  (all) 

Average  size  of  spatial  random  perturbations  .02 

Average  no.  of  timesteps  between  random  perturbations  2 

Period  of  1.68 

Mean  velocity: 

Point  1  2.20 

2  2.03 

3  2.95 

4  1.85 

Sample  size/subsample  size =  no.  of  subsamples  512/64  =  8 
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FIGURE  S15 


Bifurcation  parameter  n  .325 

No.  of  timesteps  500 

Stepsize  .03 

Init.  condition  for  0  (all) 

Average  size  of  spatial  random  perturbations  .015 

Average  no.  of  timesteps  between  random  perturbations  2.5 

Period  of  1.66 

Mean  velocity: 

Point  1  1.95 

2  1.64 

3  2.65 

4  1.60 


Sample  size/ subsample  size  =  no.  of  subsamples 


512/64  =  8 


E./U.HRSS*URVEtf  E./U.mSS*WRVEtf 
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FIGURE  S16 


Bifurcation  parameter  n  .33 

No.  of  timesteps  645 

Stepsize  .03 

Imt.  condition  for 

2  rest  =  0 

Average  size  of  spatial  random  perturbations  .01 

Average  no.  of  timesteps  between  random  perturbations  2 

Period  of  1.65 

Mean  velocity: 

Point  1  1.75 

2  1.58 

3  2.43 

4  1.59 

Sample  size/subsample  size  =  no.  of  subsamples  448/64  =  7 


(First  200  timesteps  are  unperturbed.) 
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FIGURE  S17 


Bifurcation  parameter  n  .445 

No.  of  timesteps  500 

Stepsize  .03 

Init.  condition  for  .1  (all) 

Average  size  of  spatial  random  perturbations  .01 

Average  no.  of  timesteps  between  random  perturbations  2 

Period  of  1.55 

Mean  velocity: 

Point  1  1.61 

2  1.68 

3  2.67 

4  1.76 


Sample  size/subsample  size  =  no.  of  subsamples 


512/64  =  8 


E./U.MRSS*URVEtf  E./U.MRSS*URVE# 


S17 
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APPENDIX  G 


1.  PROGRAM  LISTINGS 

The  following  programs  and  files  are  listed: 

1)  DSCRB  (+  subroutines) 

2)  ADV 

3)  BIF3A  (+  subroutines  FUN,  FUNP,  MFC) 

4)  ODE.M  (+  subroutine  F45) 

5)  Control  parameter  file  4. NS. A. 

DSCRB  with  subroutines  and  ADV  are  described  in  Appendix  B. 

BIF3A  is  the  calling  program  for  the  subroutine  BIFOR2  (see  Chapter  8). 

BIFOR 2  is  not  listed  here,*  since  it  could  be  used  in  its  original 

form  with  only  the  following  modifications  for  the  present  study: 

1)  Conversion  to  double  precision. 

This  was  done  by  W.G.  Aiello,  Department  of  Mathematics, 
University  of  Alberta.  It  was  necessary  since  BIFOR 2  was  written 
on  a  CDC  6400,  which  is  in  its  single-precision  mode  almost  as 
accurate  as  an  AMDAHL  (computer  at  U.  of  A.)  in  double  precision 
mode.  Therefore,  algorithms  accumulated  roundoff  errors  too  fast 
and  could  not  meet  their  tolerances  on  AMDAHL  single  precision  mode. 
This  was  the  only  major  modification. 

2)  The  switch  ICK  was  added  (via  parameter  string)  to  skip  the  routine 
CHECKJ  on  option,  which  evaluates  the  Jacobian  of  the  system  by 

A  listing  is  available  as  microfiche  in  [5]  or  directly  through  its 
author:  B.  Hassard,  Dept,  of  Mathematics ,  SUNY,  Buffalo,  NY. 
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finite  differencing  and  compares  it  against  an  analytic  evaluation 
done  by  the  user-supplied  subroutine  FUN  (or  FUNP  respectively) . 
This  is  just  a  CPU-time  saving  measure,  since  this  test  has  to  be 
performed  only  once  to  insure  a  correct  analytic  evaluation  of  the 
Jacobian  by  FUN. 

3)  The  parameters  EPSR  (tolerance  for  the  above  mentioned  comparison 
test  and  zero-criterion  for  the  eigenvalue  of  the  Jacobian) ,  EPS, 

NSIG  (parameters  used  for  zero-criteria  in  various  algorithms  and 
iterations  used)  and  ITMAX  (maximum  number  of  iterations  allowed 

in  various  algorithms)  can  be  set  externally  through  parameter  string. 

4)  The  microblock  COMMON/CT/ICT  was  introduced  to  BIFOR 2  and  FUN. 
The  switch  ICT  is  reset  from  zero  to  one  after  the  critical  value 
of  the  bifurcation  parameter  has  been  found  and  is  used  in  FUN  to 
save  a  copy  of  the  coefficient  matrix  of  the  system  and  its 
transpose  at  that  point  for  evaluation  of  the  mean  flow  in  MFC 
and  possibly  the  second  and  third  level  bifurcations  in  FUNP. 

A  peculiarity,  which  is  MTS*-specific ,  is  the  requirement  to 
multiply  the  line  number  in  sequential  READ  and  WRITE  statements  by 
a  factor  of  1000.  This  can  be  seen  in  the  transposing  sections  of 
BIF3A  and  ADV  and  in  all  READ  statements  from  logical  unit  4 
(attached  to  control  parameter  file  4. NS. A). 

To  set  up  a  particular  discretization  and  flow  profile,  one 
would  proceed  as  follows: 


•k 

Michigan  Terminal  System:  Operating  system  used  on  the 
University  of  Alberta,  AMDAHL  computer. 
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1)  Set  maximum  number  of  increments  in  X-,  Y-,  Z-direction  on  line  136 
of  4. NS. A.  The  numbers  chosen  are  equal  to  the  numbers  of 
corners  of  rectangular  blocks  and  equal  the  numbers  of  blocks  minus 
one  respectively  in  the  particular  directions  (maximum  number  of 
increments  is  5) . 

2)  Set  corner  coordinates  of  blocks  in  X- ,  Y-,  Z-direction  in  lines 
138  to  140  of  4. NS. A  (maximum  number  of  blocks  per  direction 
is  4)  . 

3)  Set  powers  of  the  three  spatial  dependencies  for  the  X- ,  Y-,  Z- 
component  in  lines  143  to  145  of  4. NS. A.  Their  sum  per 
component  must  not  exceed  6. 

4)  Set  selection  switches  for  the  selected  terms  in  line  2  and  the 
values  of  their  coefficients  in  line  3  of  4. NS. A.  If  affine  or 
power  dependencies  on  the  bifurcation  parameter  n  of  these 
coefficients  are  desired,  set  the  appropriate  constants  CC  and 
powers  IX  in  lines  6  and  8  respectively. 

5)  Set  the  dimensions  of  ACIR,  ACIT  in  BIF3A  and  FUN  equal  to 
the  number  of  terms  selected  in  Step  4) .  The  lines  which  require 
this  operation  are  marked  by  a  star  parade. 

6)  In  subroutine  FUN  set  up  the  final  flow  profiles  in  the  section 
"setting  of  flow  structure",  again  marked  by  a  star  parade.  Make 
sure  the  coefficients  CF,  as  selected  in  Step  4),  are  arranged 
in  ascending  order  w.r.t.  their  indices. 

The  setup  in  the  listings  presented  is  for  the  study  of  flow 

IV  (see  Chapter  9) . 


j 
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The  graphics  and  Fourier  analysis  programs  are  rather  straight¬ 
forward  and  technical  in  nature  without  any  details  of  interest  w.r.t. 
the  problem  under  study.  They  are  therefore  not  listed.  The  same  goes 
for  the  transposing  routine  TRNSP  (see  Appendix  B) ,  since  it  is 
identical  in  its  function  and  structure  to  the  transposing  section  in 
the  program  ADV. 

The  following  program  chart  (Figure  17)  shows  the  sequence  of 
execution  of  the  programs  with  logical  unit  numbers  for  data  transfer. 
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187  .  C 

188.  C  SIXTH  (INNERMOST)  LOOP  (L6)  OF  SIX  NESTED  LOOPS  . 

189.  C  ITERATING  THROUGH  THE  FOUR  TERMS  OF  BARYCENT.  COORDS. 

190.  C  EXPRESSED  IN  CARTES.  COORDS. 
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134  PARAMETERS  FOR  DISCRETIZATION  PROGRAM  DSCRB 

135  1  . DAT  A  LINE:  MAX. NO. OF  INCREMENTS  IN  X-.Y-.Z-  DIRECTION 

136  222 

137  2. -4. DATA  LINE:  DISCRETIZATION  INCREMENTS  IN  X- , Y- , Z -DIRECT  I ON 

138  000.000  001.000  010.000  003.000  004.000 


